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ABSTRACT 


Model building is a very important process in any scientific investigation. The 
researcher requires a mechanistic model, which is capable of describing the process 
completely. The final form of the mechanistic model may be obtained by following 
several steps carefully. The first and the foremost step, viz., the estimation of 
parameters is well recognized and the next is model selection. It is a process to obtain 
the unknown parameters in the model with the greatest amount of precision possible. 
One way to estimate the parameters for a large variety of models is to use the least 
squares technique, which involves minimizing the residual sum of squares of the 
differences between the experimental and the model values. In large majority of real 
systems, the investigator frequently comes across situations where several models 
have been postulated for the same system. This necessitates the discrimination among 
the competing models and selection of the best model. 

In the present investigation the Asymptotic theory and Bootstrapping methods 
are employed for the parameter estimation process in the single response as well the 
multi response systems. Initially simulated data are used for the estimation of the 
parameters and the efficiency of the above mentioned methods are tested. The model 
selection is also carried out by using Cross Validation procedure. Using a real life 
data the Asymptotic theory and Bootstrapping methods are again used for parameter 
estimation in single response models and the estimates of the parameters are 
compared with the values given by the investigators. 


XI 



INTRODUCTION 

CHAPTER 1 

INTRODUCTION 


In the study of a system one often wishes to know the behavior of the process. 
The best way to analyze a process is to put it into the mathematical equations, which 
can suitably describe the underlying phenomenon. This process is known as 
mathematical model building, which is a widely used tool and not untouched by any 
stream. In model building there is a dependent variable yi (whose true value is Tji) 
related to a set of independent variables x = (xi, X2, X3, . . .) and a vector of parameters 
p = ( pi, P2, P3...)' of the process under study such that the equation shows the true 

nature of the process. Any phenomena can be represented mathematically by 

= ( 1 - 1 ) 
where p is the observed response, x is the vector of independent variables and p is 
the vector of parameters in the equation and t] is related to x and p by ftmction f. 

This is known as the true mathematical model. In practice, the true 
mathematical model of a process is never precisely known. In many situations, the 
mechanism of a process is not completely understood and a very close mechanistic 
mathematical model representation is available. On the other hand, many situations 
occur where the researcher often wishes to obtain an approximate idea of the behavior 
of the response over some region of interest. When the mechanism is complex and its 
knowledge involves great effort, it may be wasteful to attempt to determine a 
mechanistic modely. In such cases the empirical approach of model building is used 
which does not need the knowledge of the mechanism. The response or dependent 
variable may safely be taken as a polynomial f ( x , P ) in the smooth region of interest 

y = f(x,p) = po + PlXi + P2X2+ ... + pkXk+pi2XiX2+... 

+PiiXi^+...+PkkXk^ +e (1.2) 

where p is a vector of parameters, x is the vector of independent variables, y 

is the observed response and e is the vector of error values. 
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The experiments are conducted in a region of interest and an adequate 
polynomial fitted by the least square method. It is not suggested that an empirical 
approach is used in cases where the objective is to know the mechanism, when this 
can he done simply hy mechanistic modeling. 

1.1 Uses of Mathematical Modeling 

The mathematical modeling process is generally used in each and every 
scientific investigation. It is a very useful tool to compress large amoxmt of 
information, providing new insights into the process and suggesting ways for the 
future development of the system. 

Often there exists a functional relationship, which is too complicated to grasp 
or to describe in simple terms. In this case, the experimenter may wish to approximate 
this functional relationship by some simple mathematical function, such as a 
polynomial, which contains the appropriate variables and which graduates or 
approximates to the true function over some limited ranges of variables involved. By 
examining such a graduating function, the investigator may be able to learn more 
about the underlying true relationship and to understand the effects produced by 
changes in certain important variables. 

It becomes possible by model building to predict responses and optimize the 
system. It is very important in industrial processes when it is very important to have 
scientific knowledge of what is actually happening. The design, optimization, control, 
maintenance etc. all steps become simpler by analyzing mathematical equation. 

1.2 Process of Model Building 

The process of modeling is an important operation and to be carried out 
carefully in stages. Box and Hunter [7] and Kanodia [30] also represented these steps 
by the following figure. 
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Fig 1.1 Modeling Strategy 


The work of model building can be categorized in the two areas. (1) Parameter 
estimation and (2) Model specification. 


1.2.1 Parameter estimation 

The experimenter conducts some preliminary experiments and collects some 
data from these experiments. If there happens to be a single model considered 
appropriate for describing the process of the system, the only job of the experimenter 
is to calculate the precise estimates of the parameters i.e. the unknowns appearing in 
the model equation. 

A chemical engineer deals regularly with the experimental data in studies of 
chemical reactions, mixtures and processes. If the rate-determining step is known a 
mechanistic model can be derived from this information. The rate equation involves 
several unknown parameters like reaction rate constants and adsorption equilibrium 
constants depending upon the rate equation. 
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The method of least squares (applied regression analysis) is well recognized 
for the estimation of parameters. The unknown parameters are estimated in the light 
of experimental data. Suppose we have n observations (x; , yO, i = 1, 2,...n, from a 
fixed regression model with a known fimctional relationship f Then 

yi = f(Xi.P*) + ei= (1-3) 

where it is assumed that E (si) = 0 and v(f ) =I . The least square estimates of P can 

be obtained by minimizing the error sum of squares, S{fi) given by 

S(A) = E[;',-/fe,«]"- (1.4) 

/=1 

By minimizing this residual sum of squares with respect to p with the aid of 

one of several optimization techniques, the experimenter can estimate the least 
squares estimate of the parameter p . The geometrical interpretation of theory of least 
squares estimation can be better understood by representing the equation in vector 
notation. 

y=f(X,^) + f (1.5) 

where y and f are (n x 1) vectors and X is (n x k) design matrix and f is linear 
function. 

This method minimizes with respect to /?, i.e. to minimize 

i 

S = (y-:p)''(y-^) . (1.6) 

Or 

s = l|y-:?^l|. (1.7) 

2 

y-5^ (the square of the length of y-]^) will be minimum for 

A 

XP-Xp , when (y-!^)will be orthogonal to X space. In mathematical 
representation 

X^(y-:?^) = 0 (1.8) 
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xV-x^:^ = o 

;^ = (X^X)-’XV (1.9) 

This will give the least square estimates of J3 . 

1.2.2 Model Specification 

In the previous section it was assumed that the correct form of the model was 
known. Frequently in many situations involving chemical kinetics the underlying 
physical theory suggests several rival models are capable of explaining the 
mechanism of the system. In such situations, a model can be actually put to use only if 
a single model is picked out of the given lot and fitted precisely to the available data. 
This problem is commonly referred as model-discrimination problem. 

In order to reduce a large number of models to a manageable group the 
process generally used is known as model screening. Two widely used methods are 
lack of fit and residual plots. Linear models showing significant lack of fit can be 
eliminated simply by F-test on residual sum of squares. The other method, which uses 
residual plots can be employed to indicate trends peculiar to the model. 

The methods used for discrimination criteria can be categorized as Bayesian 
methods and the classical methods. The Bayesian method is based upon the Bayes 
theorem. A value is assigned to the prior probability of each of the model and then 
expected likelihood of each model is calculated and on comparing these values the 
models are discriminated. The classical methods use likelihood ratio approach. The 
different information theoretic criteria with proper choice of penalty fimction are 
widely used to choose the correct model. 

Sometimes the experimenter can not discriminate between the two equally 
good rival models firom the existing data. In this case the experimenter is required to 
carefully design the experiments to distinguish among the rival models and the design 
of experiments should be in the region of interest, where the discrimination among the 
rival models is maximum. 

It may happen after screening and discrimination analysis that a model is 
chosen that shows specific weaknesses under the certain special conditions. To correct 
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for these weaknesses the model should be modified until a form of the model is found 
that is adequate over the operable conditions for system under study. 

1.3 Single and Multiresponse Models 

Generally the true value of the response is given by the Equation (1.1). 
Suppose y is the observed response, x is the vector of independent variables and is 

the vector of parameters in the equation and y is related to x and p by function f, 

then 

y = f(x,;5) + e, (1.10) 

where e is the error associated with the measurement. If for a particular system only 
one dependent variable is observed, the system is known as single response system. 

However, situations are possible when the experimenter observes more than 
one response in a system. As in chemical reaction engineering, kinetic investigation 
frequently involves with the models where more than one rate equation is possible for 
a reaction. In such situations the experimenter has to observe several responses 
instead of one and this multiresponse model can be represented by 

Y=f(X,^) + f, (1.11) 

where the i-th response (1 < i < r) for the u-th experiment (1 < u < n) may be 
denoted by 

y? (x..«+£f. (1.12) 

It is assumed by Rao and Iyengar [36] that the errors are such that 

E(f ® ) = 0 for all i, u 

E(5®5® ) = 0 for all i, j, U9i V (1.13) 

£(£?£?) = <7, for i,j,u 

The vector =(yu*jyu'—yu^V r responses for the u*** experiments has a 
symmetric covariance matrix X given by 

S=((^!/))- i= l,2,...r&j= l,2,...r. (1.14) 
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1.4 Objective of the present investigation 

The objective of the present investigation is to devise a reliable method for the 
estimation of imbiased parameters in both single response and multi-response models. 
Initially simulated data are used for estimation of parameters from a known model for 
both single response and multiresponse models and the differences between Point 
estimation methods and Interval estimation methods for parameter estimation has 
been tested. Then the method is extended to real life problems discussed by carr [11] 
and compare with the values of parameters obtained by the investigators. 

1.5 Thesis Organization 

This thesis consists five chapters of which the current chapter is the first. The 
next chapter, Chapter 2, presents a brief literature review related to the present study. 
The methods used in the present study are explained in Chapter 3. In Chapter 4, the 
results and discussion of present work is given. Conclusions and recommendations for 
future work are reported in Chapter 5. Finally the references used throughout the 
thesis are given in alphabetical order. 


7 



LITERATURE REVIEW 

CHAPTER 2 

LITERATURE REVIEW 


In the past few years most of the studies were made on the model 
discrimination and design of experiments than parameter estimation. But the 
parameter estimation is also very important in model discrimination and design of 
experiments. Estimation in algebraic equations, linear in parameters, is well known, 
and elementary computer packages contain all the associated statistical tests. 

Though attempts have been made to find out the reaction kinetics for different 
types of reactions, selection of objective function for nonlinear system has always 
been a matter of contention. Various methods for estimation of parameters in case of 
single and multi response systems are available in literature, but the choice of using 
objective function for nonlinear regression is quite interesting. Also there are number 
of optimization methods are available for the estimation of parameters by minimizing 
objective function. In recent years the breakthroughs in statistics is the using of 
bootstrapping methods for regression problems. Since, many real world applications 
are begiiming to appear along with numerous Monte Carlo studies on the performance 
of the bootstrap and its competitors for various problems. It is also becoming clearer 
that the bootstrap has significant in practical value but it also has some limitations. A 
summary of the literature survey on bootstrapping methods and parameter estimation 
in the area of chemical kinetics is given below. 

2.1 REVIEW OF THE PREVIOUS WORK : POINT ESTIMATIONS 

Estimation in algebraic equations, which are non-linear in parameters and in 
differential equations, was reviewed with particular emphasis on application to 
kinetics, by Seinfeld [39], Bard and Lapidus [1] and Froment and Bischoff [24]. An 
extensive treatment of the estimation problem in nonlinear regression can be found in 
the books by Draper and Smith [16a] and Bard Y [2]. Rohatgi [37] and Lehmann [32] 
have also given useful contributions in this field. Also parameter estimation in design 
of experiments for model discrimination are available in literature Mezaki et. al. [33], 
Prasad and Rao [36] etc. 


LITERATURE REVIEW 

Invariably, in several practical situations, one encounters both single and multi 
response models for parameter estimation in nonlinear regression. Seber and Wild 
[38] gave numerical methods to estimate the parameters by least square estimation, 
one of the well-recognized method is Gauss - Newton method. 

Carr[ll] studied the catalytic isomerization of n-pentane to 2-methylbutane 
and the experimental results are analyzed by applying linear least squares regression 
by transforming the nonlinear models to linear. Jhonson et. al. [29] studied the 
reexamination of hougen - watson rate models including a weighted least squares 
approach and found the dual site model is capable for explaining the data than singl 
site model in the study made by Carr[l 1]. 

Another real life data for single response models is studied by Stein [46]. He 
considered the response model of the polynomial type and parameters estimated at 
different temperatures. 

Buzzi and Forzatti [9] studied model discrimination by parameter estimation 
in five kinetic models for the synthesis of methanol from carbon monoxide and 
hydrogen. Their discrimination is based on the values of posterior probabilities for 
each model. 

Hill [26], Singh [45], studied the parameter estimation and model 
discrimination among exponential models which represents the rates of a chemical 
reaction. Kanodia [30] recently compared the efficiency of algorithms gauss- newton 
and simplex methods in parameter estimation of nonlinear models. 

Box and Drapper [6] have proposed that for multi response systems, in 
absence of variance-covariance matrix, minimization of determinant is the ideal 
objective function. They have added that the least square minimization should be used 
when the, variance-covariance matrix is known. 

Boag [5] et al. has used the minimization of the determinant as the objective 
function for a multi response nonlinear system for the oxidation of a-xylene by 
vanadia catalyst. The minimization of the determinant was achieved by using Powell's 
method. 

Vajda and Valko [48] have proposed that for multiresponse problems 
minimization of the determinant, is the objective function. But, they have observed 
that least square minimization can be used to get good initial guess or starting point 
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for obtaining the final parameters, which involves minimization of determinant. They 
have shown that least square minimization can be used in absence of correlation 
between responses; otherwise the minimization of determinant provides the best 
estimate of parameters from multi response data. 

Buzzi et. al. [10] are proposed the weighted nonlinear regression technique for 
the estimation of parameters in multi response models which is used in the sequential 
experimental design for model discrimination. 

Recently Negi [35] studied the effect of simplex method for the estimation of 
parameters in single and multi response models using nonlinear regression. 

Moros et. al. [34] described the application of genetic algorithms(GA) for 
generating the initial parameter estimations for kinetic models of catalytic processes. 
They also made a study on the rate of convergence of genetic algorithms in parameter 
estimation. 

Bates and Watts [3] fitted the Hougen-Watson model for reaction kinetics by 
using Gauss - Newton algorithm Avith Levenberg - Marquardt modifications for 
global convergence. They showed this combined procedure ensures the global 
optimum like the other methods which didn’t use differentiation for e.g. Simplex 
method, Powell's method etc. 

2.2 REVIEW OF THE PREVIOUS WORK : INTERVAL ESTIMTIONS 

Although regression analysis is one of the most widely used statistical 
techniques, application of the bootstrap to regression problems has only appeared 
fairly recently. Only a few books like Sen and Srivastava [40], Draper and Smith 
[16b] and Shao and Tu [43] incorporated a brief discussion of the bootstrap in 
regression. 

Recently Efron and Tibshirini [21] and Chemick [12] gives an updated 
description of bootstrapping methods on wide range of applications. 

Early discussion of the two methods of bootstrapping in the nonlinear 
regression model can be found in Efron[17] . Efron and Tibshirani [20] provide a 
variety of interesting applications and some insightful discussion of bootstrap 
applications in regression problems. 

Bootstrapping the residuals is an approach that also can be applied to 
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nonlinear models. Shimabukuro et. al. [44] presented an early example of a practical 
application of a nonlinear regression problem. The first major study of the bootstrap 
as applied to the problem of estimating the standard errors of the regression 
coefficients by constrained least squares with an unknown, but estimated, residual 
covariance matrix can be found in Freedman and Peters [23]. 

Bickel and Freedman [4] also deal with issues regarding bootstrapping in 
regression problems. Their study is very interesting because it shows that the 
conventional asymptotic formulas, which are correct for very large samples, do not 
work well in small to moderate sample size problems. They show that these standard 
errors can be too small by a factor of nearly three. On the other hand, the bootstrap 
method gives accurate answers. 

Theoretical work on the use of bootstrap in regression is given in Freedman 
[22], Weber [49], Wu [50] and Shao [41,42]. A subset of articles that present the 
theoretical justification for the bootstrap are Efron [18,19]. 

These authors pointed out that there is unfortunately no good rule of thumb to 
apply to determine when the conventional formulas will work or when it may be 
necessary to resort to the bootstrap. They suggest that the development of such a rule 
of thumb could be the result of additional research. Even the bootstrap procedure has 
problems in this context. 

Efron and Tibishirani [21] discuss the methods called the parametric bootstrap 
and non-parametric bootstrap. The parametric bootstrap is closer to the Monte Carlo 
methods. They also discussed about bootstrap estimate of bias. 

Breimn [8] gave various recent regression applications and also the model 
selection using bootstrapping methods with the help of Cross validation technique. 

Hayes et. al. [25] have extended bootstrap methods to the case of several 
unrelated samples with application to estimating contrasts in particle physics 
problems. 

2.3 REVIEW OF THE PREVIOUS WORK IN MODEL SELECTION 

In many situations involving chemical kinetics the underlying physical theory 
suggests several rival mechanisms that give rise to different models. The problem 
becomes one of the choosing the best model from among a group of ‘m’ rival models. 
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Cox [13,14] used an hypothesis testing approach to choose between two rival 
response functions. Cox considered the weighted sum of the functions and tested the 
null hypothesis versus the alternative hypothesis. 

Hunter and Reiner [28] suggested a sequential procedure for design of 
experiments which leads to model discrimination. They suggested the set of operating 
conditions that maximized the deviation between the predicted values of the response 
was chosen for the next nm. 

Recently, methods like AIC, BIC and EDC are proposed for the choice of a 
model by minimizing a criterion function defined on the set of alternative models. 
Hocking [27] and Thomson [47] studied the above different methods in the case linear 
regression problems. Recently, Kundu D and Murali [31] are discussed about the 
model selection in linear regression by using Monte Carlo simulation. 

2.4 SUMMARY 

From the above review, it is clear that parameter estimation and model 
selection is very important in kinetic modeling. The literature review of interval 
methods shows the importance of applying them to nonlinear models. To understand 
the effect of these methods on parameter estimation and cross validation methods for 
model selection, a simulation study is carried in this thesis. The data is simulated 
using nonlinear models, which generally observed in chemical reactions. The effect of 
parameter estimation for different variances is also studied in both single and multi 
response models. The parameters are estimated by using Gauss-Newton method with 
Levenberg-Marquardt modifications. 
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CHAPTER 3 

METHODOLOGY 


This chapter deals with the techniques and methods used in the parameter 
stimation and model discrimination for both the single response and the 
lultiresponse models. There are two ways in which one can get the data - either using 
jal life systems or by generating the data by simulation. In the present study the data 
re generated using simulation. The data simulation is done for single response as 
»rell as for multiresponse example. In addition, a real life example has been 
onsidered later. 

J.l Single Response Models 

For a particular system only one dependent variable is observed then the 
system is said to be single response system. Suppose y is the observed response, x is 
the vector of independent variables and is the vector of parameters in the equation 

and y is related to x and by function f, then a single response model is given by, 

y = f(x,;^). (3.1) 

For example, suppose a chemical reaction is occurring in the following way 

A ► B. (3.2) 

Depending on the order of the reaction, the expected concentration of A is given by, 

(3.3) 

where ti is expected concentration of A. 

01 and 02 are the parameters involved in this model. 

and ^2 are the reaction time and the temperature respectively. 

Here equation (3.3) represents a single response model for a chemical reaction 
used by Hill [26]. This example consists of a nonlinear equation of exponential type. 
There are various single response models in chemical engineering of different types 
like polynomials, differential equations and other nonlinear equations. 
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3.1.1 Simulation of the Data 

The data for and ^2 2 nd the values of parameters are taken from Hill 
W.J.[26]. Using these values the theoretical concentration of A are calculated as 
follows. 

y = {1 + <9,^. exp(-0^ / 

R‘ = R*a + p (3 5) 

where y is the generated value of concentration of A, r| is the theoretical (true) 
value of concentration of A, R* is the random number generated at a mean (p) and 
variance (cr^). R is the random number generated from a random population, which is 
distributed as N (0,1). The response data is generated at different variances by all 
values of ^1 and ^2 and the parameters are estimated. 

3.1.2 Parameter Estimation 

The method of least squares estimation is used to estimate the parameters. The 
objective function is defined as 

= + . (3.6) 

e ^,02 /=! 

Using the simulated values of concentrations, this function is minimized with 
respect to 0i and 02 with the help of Guass Newton-Levenberg Marquardt method, 
which is described in the section 3.3. The 95% confidence intervals for the parameters 
are calculated by simulations using Asymptotic theory and Bootstrapping methods. 
The average estimates, mean square errors of the parameters are also calculated. 
Another real life example is discussed by Carr [11], which involves the parameter 
estimation by least squares minimization. 

3.1.3 Model Selection 

Suppose an experimenter is investigating a system in which the underlying 
theory suggests at the outset that there are number of models which might represent 
the system. This is not new in chemical reaction engineering. For example, if a 
chemical reaction is observed as given in equation (3.2), depending on whether the 
reaction is of first, second, third or forth order, we can consider four rival models. In 
this case identifying the correct model or model selection plays important role since 

the correct model gives exact order of the reaction, which is useful to experimenters 
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For example, the four rival models for equation (3.2) are given as follows. 


=exp{-6',<"f,exp(-e“/4)} 

(3.7) 


(3.8) 


(3.9) 

r‘«={l+3«,'“^,exp(-«,"> /&)}■* 

(3.10) 


Hence the model selection plays an important role in data analysis. Here 
model selection is done with the help of Cross Validation procedure, which is 
explained in the section 3.5. The data and the rival models which are considered for 
model selection are given by Hill W.J.[26]. 

3.2 Multi Response Models 

There are some situations in which the experimenter observes more than one 
response in a system. As in chemical reaction engineering, kinetic investigation 
frequently involves with the models where more than one rate equation is possible for 
a reaction. In such situations the experimenter has to observe several responses 
instead of one and this multi response model can be represented by 

y = f(X,;^)- (3-11) 

where Y is the observed response, X is the vector of independent variables and is 

the vector of parameters and f is a nonlinear function. Here the equation for i* 
response (1 < i < r) and u* experiment (1 < u < n) is represented by, 

y?=' 7 f(X.,«. (3.12) 

For example a multi response model can be represented by the following 
reaction scheme given in thesis by Negi [35]. 
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where AjBjCjD and E denote a-Pinene (C 10 H 15 ), a- and P- Pyronene (CioHie), 
Dependene (CioHie), D. a//o-Ocimene (CioHie) and E. Dimer (C20H32) respectively 
and ki, k2, ka, k.3, k4, k4,k5 are the reaction rate constants. 

The rate equations for the model shown in fig. 3.1 are given by. 


r,=-(A,+^)C,- 2 iqC/ (3.13) 

'^K^b (3 • 14) 

^c=KCa (3.15) 

(3.16) 

r,=k,Q^+k,Cj,^-k^Q (3.17) 


which generally expressed as a multi response system. Here r represents rates of 
reactions, C represents concentrations of the species and k represents the rate 
constants, which are unknown parameters. 

3.2.1 Simulation of the Data 

The theoretical (true) rates are calculated by using the values of independent 
variables. The response values are generated by adding random numbers to these true 
values which are generated at different variance-covariances. Let a nonlinear multi 
response function is given by, 

Y = f(X,^) + f. (3.18) 

here s represents the vector of random numbers generated from N(|j., S). It is 
generated as follows. 

s = Z^R+ju (3.19) 

where ^ is a vector of standard normal random numbers, 

// is a vector representing the mean, 

is a matrix such that Z = S and S is covariance matrix. 

As described in the above procedure multivariate random errors are generated 
and added to the true values generated by the rate equations. 

3.2.2 Parameter Estimation 

The parameter estimation in multi response models is not as simple as in the 
case of single response models. Since the order of the dependent variables is different 
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for different responses and also there may be correlation between the responses. 
Hence different procedures are available for parameter estimation in multi response 
models as seen in literature. The method of weighted least squares estimation is used 
to estimate the parameters for this type of problems. The parameters are estimated by 
minimizing the objective function given by, 

0-F.,='Z(ym-y,(X.A))" Q (y,(~X.)-y,(.X.A))- ( 3 - 20 ) 

u=l 

where Q represents the reciprocal covariance matrix (2). 

Here also the minimization of objective function is carried with the help of 
Guass Newton-Levenberg Marquardt method, which is described in the section 3.3. 
The 95% confidence intervals for the parameters are calculated by simulations using 
Asymptotic theory and Bootstrapping methods. The average estimates, mean square 
errors of the parameters are also calculated. A simulated study is made on this 
estimation as discussed by Buzzi et. al. [10]. 

3.2.3 Model Selection 


It is quite common in chemical engineering practice where more than one 
chemical reaction takes place and diey are typically observed by models that include 
more than one response. Multi response models which frequently arise in 
experimental situations contains nonlinear equations and capable of describing the 
mechanism of the given system equally well. Hence therefore, they are considered to 
be rival models and the identification of correct model is a challenging problem. 

For example, four chemical kinetic models, each having two responses given 
by Buzzi et. al. [10] are. 




( 1 ) 






^2,1^1 ^2 


(1 + (1 + jXj + ^^4 jXj) 


Ti 


( 2 ) 


^,2 ^1^2 


( 2 ) 


Tl 


(3) 


^ 3 X^X 2 ^ ^ 

^’^2 - 


+ i\ + K33X,f’ 


Ti 


(4) 


*,4X1X2 


^2.4^1^2 

(l + r3,4X,+J^4.4X3)'-'^ ~(l + K3^,Xy) 


.T2 


(4) 


(3.21) 

(3.22) 

(3.23) 

(3.24) 
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Here the model selection is done with the help of Cross Validation procedure, 
which is explained in the section 3.5. The data and the rival models which are 
considered for model selection are given by Buzzi et. al. [10]. 

3.3 Gauss Newton - Levenberg Marquardt Method (GNLM) for 
Parameter Estimation 

This method was first discussed by Bates & Watts [3] and they said it is a 
powerful one in searching the optimum. This method is nothing but the combination 
of two methods Guass-Newton & Levenberg-Marquardt. In this method the iterations 
of the Gauss-Newton method are carried using the Levenberg-Marquardt method to 
modify the Hessian matrix if the search direction is to far away from the direction of 
the steepest descent. Hence the failure of the Guass-Newton method will be taken care 
in most cases where the Hessian matrix is near singular. To understand this method 
properly first Guass-Newton method and later Levenberg-Marquardt method are 
explained. 

i) GAUSS-NEWTON METHOD 

Guass-Newton gave the method as described below. 

If the functional form of the model is 

Y = f(^,e) + s (or) T) = E(y) = f(^,e), (3.25) 

where f(^,0) is non-linear, then, for a single observation in the k-th experiment, one 
may write 

yic = f(^k,e) +£k. (3.26) 


and the error sum of squares is given by, 

k=1 

Since yk and ^ are known, the sum of squares is a function of 0 only. The 0 
are found such that the sum of squares of errors is a minimum. To find the least 

A ■ 

squares estimates of 0, the function is first expanded about an initial estimate in 
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Taylor’s series and truncated after the first partial derivatives. The linear resultant 
function is, 

(3.28) 

by abbreviating, 

+ (3.29) 

Using this approximation in minimizing the sum of squares with respect to 0, the next 
approximation of 0 is given as 

( 9 ( 0 + 1 ) ^^( 1 ) ( 330 ) 

where 

- 1-1 

F(0^‘'^Y F(0^‘‘^) F(0^°^f{y-f{0^°^). (3.31) 

This provides the iterative scheme of obtaining 0. This procedure is iterated until the 
correction S becomes exceedingly small. 

Seber and Wild [38] also gave a more general approach, viz., the Newton 
method, which applies to any function satisfying appropriate regularity conditions and 
in which the residual sum of squares S (0) is expanded directly using a quadratic 
Taylor expansion. The minimization of residual sum of squares occurs, when 


^(0+1) ^^(1) (-3 32) 

where 0^'^ is the initial guess of estimates and 

5^“'^ = [/f (^<"^)]'' Ji0^°^) (3.33) 

i.e. H{0^‘‘^ ) (3 .34) 
which can also be rewritten 

[/(^"Y /(^")) + 5(^"^)] = -J(^^)) . (3.35) 

Neglecting the second derivatives matrix we obtain the “normal equations” 

and the Gauss-Newton direction 

= (3-36) 

where J(0) = - (3.37) 

60 
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rhe convergence criteria used for halting the procedure is given by 


A^y 


<r. 


(3.39) 


where /is some suitably small value (tolerance) greater than zero. 

A * 

If convergence is not achieved ^ is updated by replacing the old values by 
the new values and the process repeated. 


ii) LEVENBERG-MARQUARDT METHOD 

Levenberg Meirquardt gave the method as described below. 
Statement of the Problem: The model to be fitted to the data is 


E(y) = f(Xi,X2,. 


• • I X^ 1 , 02 1 • • • 1 0fn ) 


= f(x,0, 


(3.40) 


Where , Xj , . . . , X^^ are independent variables, , ^2 > • ■ • ■ the population 
values of m parameters, and E (y) is the expected value of the dependent variable y. 
The data points are denoted by 

-^w)- ' = '*■2 n. (3.41) 

The problem is to compute those estimates of the parameters that will 
minimize 




(3.42) 


A 

Where y , is the value of y predicted by (3 .4 1 ) at the i*'' data point. 

The model is linearized by expanding y] in Taylor series about the current 

trial values for the coefficients and retaining the linear terms only. 


A 


y ' =y 


" A “ 

* 

" A ” 

* 


A 

A 

A^i + 

5y, 

A 

A 

A6 2"^ 

dYi 

A 

1 

CD 

1 


d02_ 




A0m, 


(3.43) 
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where 


A 



The asterisk designates that the quantities are evaluated at the initial trial 
values. The linearized model is substituted into the objective function and the normal 
equations formed by setting the partial derivatives of the objective function with 
respect to each coefficient equal to zero, 

-^ = 0. j = l,2, ,M. (3.44) 

dOj 

Thus h.6 j can be found by solving 


^ My = g 

Where 


pin 


dbj 

V J 


, i= l,2...,n;j = 1,2,. ..,m. 


g = 


M 


Yi-Yi 


* ^ 


dY; 


dOj 


The resulting normal equations will be of the form 


Mi+^yDMy=/(Z-Z_)’ 


(3.45) 


(3.46) 


(3.47) 


(3.48) 


Where / is an identity matrix of order m x m and Xj is a factor that is added 

to the main diagonal of the matrix. The rules for calculating A^ are discussed in 
the original article by Marquardt are as follows. 

1. Compute 

2. Pick a modest value for Xj , say Xj = 0.001. 

A 

3. Solve the linear system of equations (3.48) to find E.0 j and calculate 

/[«; +m/ ■ 

4. If /f^y+Myl > /(^^ ) , increase Xj by a factor of 10, and go back to (3). 
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5 . If / + AO j j < f[Oj^, decrease Zj by a factor of 1 0, update the trial solution, 

A 

i.e. replace 6j by Oj +A0j and go to next time step. 

The Jacobean J, and error matrices are given by, 


dyi dy^ 


dO\ d02 


J= . 


A A 

dO, 802 


{0^-0^) 
A A * 

{ 02 - 02 ) 


A0 = \ 


Vi -Vi 


72-/2 


(y-y ) = 


{0m—0m) 


Yn-Y^ 


The vector and ^ will approach zero as convergence is achieved. If convergence 
is achieved the final coefficients are calculated from 


^ =0 + A0 , 


j-l,2,...,M. 


The convergence criteria used for halting the procedure is given by 


(3.49) 


Where y is some suitably small value (tolerance) greater than zero. 
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If convergence is not achieved 9_ is updated by replacing the old values by the ne-w 
values and the process repeated. 

In GNLM procedure as said above the modification is just the replacing of 
equation (3.36) with equation (3.48) which avoids the singularity of J J matrix. A 
flow sheet illustrating the GNLM procedure is given in Fig. 3.2. 



Fig 3.2 GiVLM Method Logic diagram 
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3.4 Asymptotic Properties of Non-Linear Least Squares Estimation 

Assume the functional form of a nonlinear model is, 

y,=m + s,. (3.50) 

where t = l,2,3,...,n. 

St are independent identically distributed errors with zero mean and variance 

0 ^ > 0 . 

Then the sum of squared errors function is given by, 

s(e)^YXy-fXe)f. (s.si) 

t=\ 


A 

The least squares estimate 6 is the vector which minimizes equation (3.51) 
and the solution of g*(^) = 0. 

then, s\d)-s\e)=ie-e)s\e) ( 3 . 52 ) 


^ is a point on the line joining 6 and 6 . 


where = 


dS{e) dS{9) 


V 


d9, 


p J 


d^S(9) d^S(9) 
69^ "'d9,d9p 


S\9) = 


d^S{9) 

d9pd9^ 


andiS’(^) = 0. 


d^S(9) 

d9l 

-S\9) = {9~9)S\9) 
S\9) = -2[F{9)F{9f) 


(3.53) 

(3.54) 


We assume that [F{9)Fi9)^'\ is a continuous function of 9 and the inverse 

A 

exists near ^ or ^ . 


S-W =-^SW = -2[ WCT'] 

S \ff) = -22(3', - m)^f.(e) = -2eF(,ef (3.55) 
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The expansion is valid V ^ in a neighborhood of G and in particular for 6 
also for large n. 

-sFiGf =(e-e)[F{e)F{ef'] (3.56) 

0-9) = -£F0f (3.57) 

=^yf^0-G) = -y[^ l£F0Y [F0)F0f J 

= (3.58) 

as n is large then (o^GfuG^ and 


i[F(^)F(^)^]] ->[i[F(^)F(^)"] 


->S 


where E is a positive definite matrix. 


where, 


I i 

Np(0,CT^I) E’* 

4^0-0,) = -{X.C) ~ N, {0,C^D(X)C) 

4n0-e,) = Np (0,o-'E-'EE-‘ ) = Np [0,cr^ir' ) 

S = limi(F(«.)m)’') 


(3.59) 


(3.60) 


(3.61) 

(3.62) 


This asymptotic distribution provides us to construct confidence intervals and 
testing problems. 

The necessary assumptions which are required to prove the consistency are, 

1 . E(£i) = 0 and V(si) = 0 ^ , Si ‘s are i.i.d. 

2. f(6) is a continuous function of ^ . 

/\ A 

3. ^ € © , 0 is compact set (closed & bounded). ^ is an interior point of 0 . 

4. f(d)is continuous and bounded for all 0 . 

5 . Implies that - [/(^i ) - /(^ 2 )] converges uniformly to 00^ , 02 ). 

6 . — ^[/(^i)-/(^2)J convergesto £>(^1,^2) 
iff <9 i=^2- 
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3.5 Bootstrapping Methods for Parameter Estimation 

The treatment of the bootstrap methods described here comes from Efron and 
Tibshirani [21]. The interested reader is referred to that text for more information on 
the underlying theory behind the bootstrap. There does not seem to be a consistent 
terminology in the literature for what techniques are considered bootstrap methods. 
Some refer to the resampling techniques of the previous section as bootstrap methods. 
Here, we use bootstrap to refer to Monte Carlo simulations that treat the original 
sample as the pseudo-population or as an estimate of the population. Thus, in the steps 
where we randomly sample from the pseudo-population, we now resample from the 
original sample. 

3.5.1 General Bootstrap Methodology 

The bootstrap is a method of Monte Carlo simulation where no parametric 
assumptions are made about the underlying population that generated the random 
sample. Instead, we use the sample as an estimate of the population. This estimate is 

A 

called the empirical distribution F where each Xt has probability mass 1/n . Thus, 
each JIG has the same likelihood of being selected in a new sample taken from F . 

When we use F as our pseudo-population, then we resample with replacement from 
the original sample x = (Xi,...,jc„). We denote the new sample obtained in this manner 

by X = (x*,..., jc*) . Since we are sampling with replacement from the original sample, 

there is a possibility that some points Xt will appear more than once in x* or maybe 
not at all. We are looking at the univariate situation, but the bootstrap concepts can 
also be applied in the tZ-dimensional case. 

A small example serves to illustrate these ideas. Let's say that our random 
sample consists of the four numbers x = (5, 8, 3, 2). The following are possible 
samples x*, when we sample with replacement from x : 

= (X4.X4,X2,Xi) = C1,2A5) 

X*" = (X4,X2,X3,X4) = (2,8,3,2). 

We use the notation X*’’, b = 1,..., B for the b-th bootstrap data set. 


26 



METHODOLOGY 


In many situations, the analyst is interested in estimating some parameter e by 
calculating a statistic from the random sample. We denote this estimate by 

e = T = tix„...,x„). (3.63) 

We might also like to determine the standard error in the estimate 6 and the 
bias. The bootstrap method can provide an estimate of this when analytical methods 

A 

fail. The method is also suitable for situations when the estimator 0 = t(x) is 
complicated. 

To get estimates of bias or standard error of a statistic, we obtain B bootstrap 
samples by sampling with replacement from the original sample. For every bootstrap 

A 

sample, we calculate the same statistic to obtain the bootstrap replications of ^ , as 
follows 

b = (3.64) 

These B bootstrap replicates provide us with an estimate of the distribution of 

A 

6. This is similar to what we did in the previous section, except that we are not 
making any assumptions about the distribution for the original sample. Once we have 
the bootstrap replicates in equation (3.64), we can use them to understand the 
distribution of the estimate. 
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The steps for the basic bootstrap methodology are given here, with detailed 

procedures for finding specific characteristics of 0 provided later. The issue of how 
large to make B is addressed with each application of the bootstrap. 

PROCEDURE FOR BASIC BOOTSTRAP 

1 . Given a random sample, x = (x, , . . ., ) , calculate 0 . 

2. Sample with replacement from the original sample to get x* = (x*,...,x*) . 

3. Calculate the same statistic using the bootstrap sample in step 2 to get, . 

4. Repeat steps 2 through 3, B times. 

A 

5. Use this estimate of the distribution of 0 (i.e., the bootstrap replicates) to obtain the 
desired characteristic (e.g., standard error, bias or confidence interval). 

3.5.2 Bootstrap Percentile Confidence Interval 

An improved bootstrap confidence interval is based on the quantiles of the 
distribution of the bootstrap replicates. This technique has the benefit of being more 
stable than the bootstrap-t, and it also enjoys better theoretical coverage properties 
(Efron and Tibshirani, [21]). The bootstrap percentile confidence interval is 

(3.65) 

where 0*^"^^^^ is the a/2 quantile in the bootstrap distribution of 0* .For 
example, if a/2 = 0.025 and B = 1000, then is the 0*’’ in the 25* position of 


the ordered bootstrap replicates. Similarly, 0*^^'^’’^^ is the replicate in position 975. As 

discussed previously, some other suitable estimate for the quantile can be used. 

The procedure is the same as the general bootstrap method, making it easy to 
understand and to implement. We outline the steps below. 

PROCEDURE FOR BOOTSTRAP PERCENTILE INTERVAL 

1 . Given a random sample, x = (Xj , . . ., x„ ) , calculate 0. 

2. Sample with replacement from the original sample to get x* = (x*,...,x*) . 


3. Calculate the same statistic using the bootstrap sample in step 2 to get the bootstrap 

Aj|j » 

replicates, 0 . 

4. Repeat steps 2 through 3, 5 times, where 5 >1000. 
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5. Order the §*'’ from smallest to largest. 

6 . Calculate B. a/2 and B. (1- a/2). 

7. The lower endpoint of the interval is given by the bootstrap replicate that is in the 

B. a/2*'' position of the ordered 0** , and the upper endpoint is given by the bootstrap 
replicate that is in the B . (1 - a/2)*'’ position of the same ordered list. Alternatively, 
using quantile notation, the lower endpoint is the estimated quantile the 

upper endpoint is the estimated quantile where the estimates are taken from the 
bootstrap replicates. 

3.6 Cross Validation Method for Model Selection 

Often, one of the jobs of a statistician or engineer is to create models using 
sample data, usually for the purpose of making predictions. We must then decide what 
model best describes the relationship between the variables and estimate its accuracy. 

Unfortunately, in many cases the naive researcher will build a model based on 
the data set and then use that same data to assess the performance of the model. The 
problem with this is that the model is being evaluated or tested with data it has already 
seen. We introduce cross-validation in regression applications, where we are 
interested in estimating the expected prediction error. 

Say we have a set of data, (^i, Tjj), where denotes a predictor variable and 77 , • 
represents the corresponding response variable of size ‘n’. We are interested in 
modeling the data. Assume there are two competitive models for the available data, 
which are given by, 

=exp{-^»J exp(-e.'Vi.)}. (3.66) 

(3.67) 

The cross validation procedure applied to the above models to select the 
correct model is explained as follows. 

Assume a model and take the given set of data. Remove the i*** observation 

A 

from the data. Then we have (n-1) observations. Estimate the 0 from these (n-1) 
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fjf^ = exp{-^®^, exp(-4 “’ /^,)} (3.68) 

Repeat the above procedure for i from 1 to n and add all the squared errors as 
follows. 

i;(>7,-^,®)'=cr(l). (3.69) 

(=1 

This is Cross Validatory error due to model 1 . Repeat the same process for 
model 2. Then the model which gives minimum value of CV(k) is considered as the 
best model among the all competing models. If the Cross Validatory error values are 
nearly equal for different models then further experiments should be done which will 
discriminate the competing models. This procedure is independent on sample size (n). 
PROCEDURE FOR CROSS VALIDATION METHOD 

1 . Let a Model(k) and sample of size(n). 

2. Skip i* Observation from the sample. 

3. Estimate parameters using (n-1) Observations. 

4. Estimate squared error value for i th Observation using new parameters. 

5. Repeat the process for i = Ito n and add all the squared errors as CV(k). 

6. The model which gives minimum value of CV(k) is the correct model. 

3.7 Estimation of Average Estimates and M.S.E. of the parameters 

The average estimates and the mean square errors of the parameters are 
calculated for the single response and multi response function. The data are simulated 
in the same fashion as it is in 3.1.1 for single response and in 3.2.1 for multi response 
example. The work of simulation is repeated 1000 times and the parameters are 
estimated 1000 times. The average estimates are obtained as follows. 

Avg. est. of k^=\l (k^ f (3 .70) 

M 

The mean square estimate of the parameters is estimated as follows. 

M.S.E. of k, = -ktf (3.71) 
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CHAPTER 4 

RESULTS AND DISCUSSION 


The methods described in the last chapter have been applied and the results 
obtained are given in this chapter. 

4.1 Single Response System 

The example given in the Section 3.1 is taken and the parameters are estimated 
for this single response system. The model is 

4.1.1 Simulation of the data 

The data are simulated for the single response system. The values of 
concentrations are taken from Hill [26] and the response values are calculated using 
the equation 

7 = {1 + 400^, exp(-5000/|J}‘' . (4.1) 

The response values are tabulated as follows. 


TABLE 4.1: Single Response data 


Run 

Input Variables 

Response 

( 7 ) 

(N) 

(^i) 

(#2) 

1 

25 

575 

0.3741 

2 

25 

475 

0.7885 

3 

125 

475 

0.4271 


125 

575 

0.1067 

5 

125 

600 

0.0768 

6 

125 

600 

0.0768 

7 

50 

450 

0.7698 

8 

100 

600 

0.0942 

9 

75 

600 

0.1217 

10 

150 

550 

0.1288 

11 

■ 25 

525 

0.5777 

1 12 

150 

525 

0.1856 1 
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The errors were generated with different values of the variances and then 
added to the predicted responses shown in Table 4.1. The generated values of the 
responses so obtained from Equation (3.4) are shown in Table 4.2. 


Table 4.2 The simulated response data with different values of variances 


CT = 0 

a = 0.001 

a = 0.01 

a = 0.05 

a = 0.1 

a = 1 

0.3741 

0.3741 

0.3525 

0.3985 

0.2349 

-1.4639 

0.7885 

0.7895 

0.7834 

0.8558 

0.7243 

0.8989 

0.4272 

0.4264 

0.4151 

0.4285 

0.4432 

0.0216 

0.1068 

0.1067 

0.1070 

0.1490 

0.0726 

-1.3660 

0.0768 

0.0779 

0.0820 

0.2079 

0.0573 

0.7915 

0.0768 

0.0787 

0.0748 

0.0573 

0.0372 

0.0673 

0.7699 

0.7717 

0.7830 

0.7668 

0.8507 

1.4772 

0.0942 

0.0937 

0.1137 

0.1088 

-0.0045 

-0.8076 

0.1218 

0.1222 

0.1318 

0.0835 

0.1458 

0.1197 

0.1288 

0.1292 

0.1492 

0.1334 

0.0660 

0.2098 

0.5777 

0.5774 

0.5659 

0.5190 

0.6512 

-0.4082 

1 0.1857 

0.1857 

0.1822 

0.2363 

0.2655 

-0.2455 1 


4.1.2 Parameter Estimation 

Using the data in Table 4.1 and the response values generated as shown in 
Table 4.2 the parameters were estimated using the GNLM method subroutine, given 
in Appendix 1. For different initial guesses and different variances the estimated 
parameters are tabulated as follows. The three different initial guesses taken were 
(410,4900); (400,5000); (390,5100). If we observe the results from Table 4.3 for the 
three different initial guesses the estimated parameters are same. Hence this assures 
that we reached the global minima. Also it is clear that the estimated parameters were 
close to the true values for lower variances. 
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Table 4.3 Estimated values of the parameters 


1 For or = 0.001 * I 

I 


..l2 

SOS 

With Initial Guess(1) 

399.67 

5000.26 

9.5547E-4 

With Initial Guess(2) 

399.67 

5000.26 

9.5547E-4 

I With Initial Guess(3) 

399.67 

5000.26 

9.5547E-4 

For a = 0.01 | 


__ 

^2 

SOS 

With Initial Guess(1) 

403.06 

4997.65 

. 

0.0018 

With Initial Guess(2) 

403.06 

4997.65 

0.0018 

j With Initial Guess(3) 

403.06 

4997.65 

0.0018 1 

For a = 0.0 

5 

I 

. ^1 . 

....^2. ... 

SOS 1 

With Initial Guess(1) 

384.24 

5010.62 

0.0309 

I With Initial Guess(2) 

384.24 

5010.62 

0.0309 

I With Initial Guess(3) 

384.24 

5010.62 

0.0309 

I Forcy = 0.- 

1 

I 

. .^1 

^2 . 

SOS 1 

j With Initial Guess(l) 

415.54 

5013.79 

0.0587 

With Initial Guess(2) 

415.54 

5013.79 

0.0587 

1 With Initial Guess(3) 

415.54 

5013.79 

0.0587 1 

1 Fora =1.1 

n 

1 

. 

^2 _ 

SOS 1 

With Initial Guess(1) 

487.52 

4699.04 

7.8705 

With Initial Guess(2) 

487.52 

4699.04 

7.8705 

I With Initial Guess(3) 

487.52 

4699.04 

7.8705 


True Values: =400,^2 = 5000. 

4.1.3 Estimation of Average Estimates and MSE of the Parameters 

The objective fimction given by Equation (3.6) is minimized by the above 
method and repeated for 1000 times. At each time the errors were generated with the 
same variance and average parameter estimations are calculated. The average 
estimates and the mean square errors of the parameters were calculated using the 
Equation (3.70) and (3.71) given in the previous chapter. The values of the average 
estimates and the mean square errors of the parameters were estimated by taking the 
different values of variances. The estimated parameters for different values of the 
variances are given in Table 4.4. However, only for lower values of variances, the 
estimated parameters were close to the true values. The mean lengths of confidence 
intervals for the parameters were calculated by using asymptotic theory 
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and bootstrapping methods given in Table 4.4. The coverage probabilities are also 
given. 


Table 4.4 Simulation Results 




0.001 

0.01 

0.05 

0.1 

- 1 

Average 

1 Estimates 


399.98 

399.91 

399.54 

401.67 

429.21 


5000.02 

5000.03 

4999.85 

4988.26 

4912.66 

Mean 

Squared 

Errors 


0.1747 

17.4880 

394.99 

1.0561E4 

8.4095E4 

^2 

0.1069 

10.5472 

226.51 

0.683 1E4 

9.32 12E4 

Mean length 
of asymptotic 
confidence 
intervals 


32.9570 

332.67 

1602.92 

8615.91 

5.9432E4 

^2 

42.4182 

428.29 

2074.82 

1.0898E4 

9.4086E4 

Mean length 

1 of bootstrap 
confidence 
intervals 


1.5229 

15.2737 

72.5226 

346.56 

776.56 

^2 

1.1711 

11.7581 

55.0437 

279.41 

884.78 

Coverage 
probabilities 
lof asymptotic 
confidence 
intervals 


994 

996 

995 

997 

996 

^2 

992 

993 

996 

994 

992 

Coverage 
probabilities 
of bootstrap 
confidence 

1 intervals 


902 

905 

899 

896 

897 

^2 

904 

905 

895 

897 

891 


The average estimates of the parameters gives the more robust values of the 
estimates, i.e. the values of average estimates do not get change easily by including 
some more data points. The mean square errors of the parameters give the deviation in 
the values of estimates with the true values of the parameters. The values of average 
estimates are again very close to the tme values, when the error included in the 
response decreases. Similarly the values of mean square errors of the parameters are 
very low when small values of variances used for generating errors. Also the 
confidence intervals calculated from both die methods are in well agreement. The 
histograms and the corresponding normality plots of the estimated parameters are 
given in figures 4.1 & 4.2. From the graphs we can clearly observe that the parameters 
are normally distributed and the assumptions of Asymptotic theory are satisfied. 
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.1.4 Model Selection 

To see whether the model, from which we generated the data turns out to be 
rue, the task of model discrimination has been performed. The cross validation 
nethod is used for model discrimination. The estimates of parameters obtained by 
)oth the methods of parameter estimation have been used for model selection. In the 
iross validation method the model, having the smallest value of CV is considered as 
he best one. The value of the cross validation index was calculated using the 
Equation 3.68 and 3.69. The results are shown in the following Table 4.5. 

Table 4.5 Cross Validation Results 


a 





1 0.001 

0.0907 

1.1485E-5 

0.0594 

0.1686 

0 

0 

0.1008 

0.0022 


0.1589 

0.05 

0.1713 

0.0355 

■iKiraari 

0.1631 

0.1_ 

o.un 



0.3049 

1 1-0 

7.5511 

6.3935 

8.9297 

9.2952 


From the above results it is evident that for different values of a, the cross 
validation index for model 2 is the least among all the rival models. This shows that 
model 2 is the best among the rival models. Also the studies made by Kanodia [30] 
says that model 2 is best based on classical methods like Bayesian and the Distance 
methods. 


4.2 Multi Response System 

Buzzi et. al. [10] considered the following chemical kinetic model which is a 
dual response system. 


Ti 




(1 + + ^^4X2) 


^2 = 




(1 + K-^X^ + K^X2) 


(4.2) 


4.2.1 Simulation of the data 

The data were simulated in the same fashion as it was done in single response 
system. The data were generated from the above model with parameters = 0.1, K 2 = 


0.01, JC 3 = 0.1, K 4 = 0.01 within operability region 5.0 <x, <55.0, 5.0 < ^2 <55.0. 
The response values are tabulated as follows. 
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TABLE 4.6: Multi Response data 


Run 

Input Variables 

Response”~~~ j 

(N) 

(^i) 

(^ 2 ) 

(T.) 


1 

20.0 

20.0 

12.5000 

HaSTiTil 

2 

30.0 




3 

20.0 

30.0 

18.1818 

■MM 

4 

30.0 

30.0 

20.9302 


5 

25.0 

25.0 

16.6667 


6 

25.0 

15.0 

10.2740 


7 

25.0 

35.0 

22.7273 


8 

15.0 

25.0 


1.3636 

9 


25.0 

18.4211 

1.8421 

10 


32.8 

26.4206 

2.6421 

11 

.. . I 

55.0 


42.9078 

4.2908 

12 

10.6 


22.3372 

2.2337 

13 

16.0 

55.0 

27.9365 

2.7937 

14 

5.0 

5.0 

1.6129 

0.1613 

15 

5.0 

55.0 



16 

55.0 



0.4198 1 

17 

55.0 

55.0 

42.9078 

4.2908 

18 

10.6 

55.0 

22.3372 

2.2337 

19 

16.0 

55.0 

27.9365 

2.7937 


The data is generated by adding the errors, which are obtained by using 
different covariance matrices. The simulated data for multi response system are 
shown in Appendix I. 


4.2.2 Parameter Estimation 

Using the values of the concentrations and the simulated rates the objective 
function given by the Equation (3.20) is minimized with respect to the K\, K 2 , K 3 and 
K 4 . The subroutine for GNLM method given in the Appendix II is used for 
minimization of the objective function with initial guess (1,0. 1,1, 0.1). The parameters 
were estimated by taking different values of the variances in simulating the data. The 
results are tabulated as follows. For lower variances the parameters are near to true 
yalues. 
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TABLE 4.7: Estimated values of the parameters 


True values: Ki = 0.1, Kj = 0.01, K 3 = 0. 


Parameter Ki 


02^ 

0.5 

0.75 

1.0 

2.0 

0.5 

0.1036 

0.0861 

0.1016 

0.1038 

0.75 

0.1050 

0.0838 

0.1065 

0.1078 







and K 4 = 0 . 01 . 


Parameter K 2 


Oi , 02 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0112 

0.009 

0.0089 

0.0122 

0.75 

0.0113 

0.0085 

0.0093 


1.0 

0.0101 

0.0128 

0.0082 

0.0144 


Parameter K 3 


Bag 









■liiiiaM 

0.75 

0.1049 

0.0807 

0.1086 

0.1056 

1.0 

0.0966 

0.1145 

0.0878 

0.1156 


Parameter K 4 


.- 2^2 
CTi , G2 


0.75 

1.0 

2.0 



0.0055 

0.0097 

0.0106 

0.75 


0.0066 

0.0117 

0.0134 

1.0 

gjfiiiiaM 

0.0160 

0.0083 

0.0125 


4.2.3 Estimation of Average Estimates and MSE of the Parameters 

The objective function given by Equation (3.20) is minimized by the above 
method and repeated for 1000 times. Each time the error will be generated with the 
same variance and average parameter estimations are calculated. The average 
estimates and the mean square errors of the parameters were calculated using the 
Equation (3.70) and (3.71) given in the previous chapter. The values of the average 
estimates and the mean square errors of the parameters were estimated by taking the 
different values of variances. The mean values of estimated parameters for different 
values of the variances are given in Table 4.8 and their mean squared errors are given 
in Table 4.9. However, only for lower values of variances, the estimated parameters 
were close to the true values. The mean lengths of confidence intervals for the 
parameters are calculated by using asymptotic theory and bootstrapping methods 
given in Tables 4.10 and 4.1 1. The coverage probabilities from asymptotic theory and 
bootstrapping methods are also given in Tables 4.12 and 4.13. The histograms and the 
corresponding normality plots of the estimated parameters are given figures from 4.3 
to 4.14. From the graphs we can clearly observe that the parameters are normally 
distributed and the assumptions of Asymptotic theory are satisfied. Also the 
parameters are unbiased for lower variance and higher sample size. 
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Table 4.9 Estimates of Mean Squared Errors 


ilts for Sample Size: 10 


Parameter Ki 


2 " 

0.5 

0.75 

1.0 

2.0 


2.4621 

2.8654 

2.5642 



3.0257 

9.4737 

9.5911 

7.1177 


5.3708 1 

9.995 

1.1514 

19.3441 


Parameter K 2 


2 2 

Ol , G 2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0273 

0.0746 

0.0827 

0.0792 

0.75 

0.0372 

0.1028 

0.1061 

0.0773 

1.0 

0.0597 

0.1167 

0.0153 

0.2171 


Parameter K 3 


2 

1 

0.5 

0.75 

1.0 

2.0 


2.9116 

2.1156 

2.5124 

3.1956 


3.1823 

10.3357 

10.4623 

8.0099 


5.5289 

2.3611 

1.1676 

21.1755 


Parameter K 4 


2 2 

oi ,02 

0.5 

0.75 

1.0 

2.0 

0.5 

0.4197 

1.3549 

1.3958 

1.0619 

0.75 

0.5937 

1.5738 

1.5941 

1.0301 

1.0 

1.1008 

1.4067 

0.2631 

3.2342 


suits for Sample Size: 15 


Parameter Ki 


2 

2 

0.5 

0.75 

1.0 

2.0 


1.9311E-4 

1.8435E-4 

2.095E-4 

2.1058E-4 

» 

) , : 

3.7151E-4 

3.2982E-4 

3.5005E-4 

3.7716E-4 


4.3323E-4 

4.3682E-4 

4.6771E-4 

4.565E-4 


Parameter K 2 


2 2 

CTl , CT 2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.285E-5 

0.35E-5 

0.46 lE-5 

0.581E-5 

0.75 

0.468E-5 

0.466E-5 

0.52E-5 

0.74E-5 

1.0 

0.523E-5 

0.544E-5 

0.695E-5 

0.846E-5 


Parameter K 3 


f 2 ^ 

0.5 

0.75 

1.0 

2.0 

! 

2.4436E-4 

2.258E-4 

2.345E-4 

2.8132E-4 

5 

4.7141E-5 

4.0853E-4 

4.3512E-4 

4.7867E-4 


5.4671E-4 

5.4527E-4 

5.8396E-4 

5.7688E-4 


Parameter Ki 


2 2 

ai ,^2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.1427E-4 

0.152E-4 

0.1651E-4 

0.1832E-4 

0.75 

0.2782E-4 

0.2568E-4 

0.2618E-4 

0.2823E-4 

1.0 

0.3205E-4 

0.3216E-4 

0.3633E-4 

0.3346E-4 


esults for Sample Size: 19 


Parameter Ki 


rr 2 
0^2 

0.5 

0.75 

1.0 

2.0 

5 

1.5612E-4 

1.4653E-4 

1.685E-4 

1.7831E-4 

fS 

2.6082E-4 

2.41E-4 

2.7996E-4 

2.6434E-4 

0 

3.5555E-4 

3.5793E-4 

3.4903E-4 

3.557E-4 


Parameter K 2 


2 ^ 2 
ai ,cr 2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.206E-5 

0.254E-5 

0.261E-5 

0.28E-5 

0.75 

0.312E-5 

0.31E-5 

0.379E-5 

0.484E-5 

1.0 

0.412E-5 

0.446E-5 

0.457E-5 

0.569E-5 


Parameter K 3 


. 02 ^ 

0.5 

0.75 

1.0 

2.0 

5 

1.8509E-4 

1.904E-4 

2.0584E-4' 

2.1163E-4 

75 

3.0875E-4 

2.8455E-4 

3.3408E-4 

3.129E-4 

.0 

4.2172E-4 

4.2341E-4 

4.126E-4 

4.2138E-4 


Parameter Kt 


2 2 

CTl ,02 

0.5 

0.75 

1.0 

2.0 

0.5 

0.1227E-4 

0.1347E-4 

0.128E-4 

0.1537E-4 

0.75 

0.202E-4 

0.1873E-4 

0.2093E-4 

0.2047E-4 

1.0 

0.2702E-4 

0.261E-4 

0.2695E-4 

0.2803E-4 
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Table 4.10 Estimates of Mean Length of Asymptotic Confidence Intervals 


ults for Sample Size: 10 


Parameter Ki 


^ 2 
G2 

0.5 

0.75 

1.0 

2.0 

5 

0.2066 

0.1647 

0.1767 

0.2248 

rs 

38.4261 

129.4084 

123.1161 

100.2374 

1 

0 

74.9554 1 

71.3931 

47.1446 

383.2195 


Parameter K2 


2 2 
,CT2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0223 

0.0178 

0.0193 

0.0247 

0.75 

4.2569 

13.4259 

12.9731 

10.4108 

1.0 

8.2075 

8.0480 

4.7599 

40.1699 


Parameter K3 


C2^ 

0.5 

0.75 

1.0 

2.0 

5 

0.218 

0.1764 

0.1891 

0.2398 

75 

39.4164 

134.0697 

128.2142 

107.0142 

.0 

86.1908 

82.2281 

48.0341 

397.5716 


Parameter K4 


2 2 

ai ,02 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0879 

0.0704 

0.0758 

0.0974 

0.75 

17.0155 

52.9719 

48.6026 

37.1306 

1.0 

23.5093 

21.8926 

20.2096 

157.0165 


:sults for Sample Size: 15 


Parameter Ki 


,CT2^ 

0.5 

0.75 

1.0 

2.0 

>.5 

0.0553 

0.0615 

0.0665 

0.0888 

.75 

0.0645 

0.0695 

0.0755 

0.0981 

l.O 

0.0697 

0.0739 

0.0819 

0.1021 


Parameter K3 


i '2 

9 G 2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0619 

0.0688 

0.0745 

0.0994 

>.75 

0.0722 

0.0777 

0.0845 

0.1098 

1.0 

0.078 

0.0828 

0.0917 

0.1143 


Parameter Kz 


2 „ 2 

CTi , 02 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0066 

0.0073 

0.0079 

0.0107 

0.75 

0.0076 

0.0082 

0.009 

0.0116 

1.0 

0.0082 

0.0087 

0.0097 

0.0122 


Parameter K4 


2 _ 2 

Oi , O2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0153 

0.0171 

0.0185 

0.0247 

0.75 

0.0179 

0.0192 

0.0209 

0.0272 

1.0 

0.0193 

0.0205 

0.0227 

0.0283 


iesults for Sample Size: 19 


Parameter Ki 


2 ^ 2 

l 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0488 

0.055 

0.0589 

0.077 

0.75 

0.056 

0.0607 

0.0669 

0.0845 

1.0 

0.0618 

0.0650 

0.0702 

0.0890 


Parameter K3 


2^2 

1 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0532 

0.0599 

0.0642 

0.0839 

0.75 

0.061 

0.0662 

0.0729 

0.092 

1.0 

0.0675 

0.0709 

0.0765 

0.0970 


Parameter K2 


^ 2 „ 2 

CTi , O2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0057 

0.0064 

0.0069 

0.009 

0.75 

0.0065 

0.007 

0.0077 

0.0098 

1.0 

0.0072 

0.0075 

0.0082 

0.0103 


Parameter K4 


2 2 

ai ,ct2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0135 

0.0152 

0.0163 

0.0213 

0.75 

0.0155 

0.0168 

0.0185 

0.0233 

1.0 

0.0170 

0.0180 

0.0194 

0.0246 
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Table 4.11 Estimates of Mean Length of Bootstrap Confidence Intervals 


ilts for Sample Size: 10 


Parameter Ki 



0.5 

0.75 

1.0 

2.0 

) 

4.9769 

4.6643 

4.5620 

3.9466 

5 

9.1990 ' 

11.11 

10.6483 

12.0243 

) 

15.3921 

16.1539 

15.2691 

16.2877 


Parameter Kj 


ai , (72 

0.5 

0.75 

1.0 

2.0 

0.5 

0.4933 

0.4651 

0.4503 

0.3855 

0.75 

0.9027 

1.1088 

1.0436 

1.1666 

1.0 

1.5293 

1.6241 

1.5507 

1.6510 


Parameter K3 


C72^ 

0.5 

0.75 

1.0 

2.0 

5 

5.5242 

5.1866 

5.0869 

4.3749 

^5 

10.1823 

12.3033 

11.7022 

13.1788 

0 

16.9915 

17.7304 

16.6584 

17.9904 


Parameter K4 


2 2 

CTi ,(72 

0.5 

0.75 

1.0 

2.0 

0.5 

1.7845 

1.6541 

1.5757 

1.3743 

0.75 

3.3336 

4.1029 

4.0012 

4.5318 

1.0 

5.9110 

6.0993 

5.8796 

6.2680 


suits for Sample Size: 15 


Parameter Ki 


,C72^ 

0.5 

0.75 

1.0 

2.0 

1.5 

0.1382 

0.0851 

0.1333 

0.0937 

.75 

0.319 

0.4481 

0.4238 

0.3457 

.0 

0.6211 

0.9057 

0.9622 

0.7749 


Parameter K2 


2 2 

ai ,(72 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0145 

0.0096 

0.0153 

0.0114 

0.75 

0.0319 

0.0442 

0.041 

0.0357 

1.0 

0.0629 

0.0936 

0.0945 

0.0731 


Parameter K3 


»cr2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.1640 

0.1003 

0.1596 

0.1109 

1.75 

0.3818 

0.544 

0.5043 

0.4147 

1.0 

0.7583 

1.0879 

1.1732 

0.9507 


Parameter K4 


2 2 

C7l ,02 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0362 

0.0234 

0.0356 

0.0256 

0.75 

0.083 

0.1137 

0.1103 

0.0898 

1.0 

0.1566 

0.2273 

0.2359 

0.1871 


iesults for Sample Size: 19 


Parameter Ki 


2 rr 2 
1 » O2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0580 

0.0585 

0.0568 

0.077 

0.75 

0.08 

0.0756 

0.1432 

0.0811 

1.0 

0.1525 

0.0957 

0.2342 

0.2613 


Parameter K3 


r 2 _ 2 

>02 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0658 

0.0666 

0.0647 

0.0839 

0.75 

0.0912 

0.0857 

0.1657 

0.0925 

1.0 

0.1782 

0.1094 

0.2641 

0.2472 


Parameter K2 


^2^2 
ai ,02 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0065 

0.0069 

0.0070 

0.009 

0.75 

0.0086 

0.0084 

0.0136 

0.0103 

1.0 

0.0164 

0.0103 

0.0257 

0.0233 


Parameter K4 


^ 2 „ 2 
» 02 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0155 

0.0155 

0.0151 

0.0213 

0.75 

0.0211 

0.0201 

0.036 

0.0213 

1.0 

0.0405 

0.0253 

0.0641 

0.0554 
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Table 4.8 Mean Parameter Estimates 


le values: K1 = 0.1 K2=0.01 

lults for Sample Size: 10 


Parameter Ki 



0.5 

0.75 

1.0 

2.0 

5 

0.1082 

0.1131 

0.1132 

0.1135 

15 

0.1599 

0.3086 

0.311 

0.3137 

0 

0.2649 

0.2603 

0.4373 

0.5141 


Parameter K3 



0.5 

0.75 

1.0 

2.0 

.5 

0.1085 

0.1137 

0.1139 

0.1142 

75 

0.1668 

0.3181 

0.3206 

0.3234 

.0 

0.2877 

0.2822 

0.4535 

0.533 


K3 = 0.1 K4 = 0.01 


Parameter K2 


2 2 
csi , cr 2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0108 

0.0114 

0.0114 

0.0115 

0.75 

0.0163 

0.0313 

0.0317 

0.0323 

1.0 

0.0279 

0.0279 

0.0449 

0.0541 


Parameter K4 


CJi ,CS2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0135 

0.0153 

0.0153 

0.0154 

0.75 

0.0309 

0.094 

0.0949 

0.0961 

1.0 

0.0635 

0.0618 

0.1454 

0.1757 


:sults for Sample Size: 15 


Parameter Ki 



0.5 

0.75 

1.0 

2.0 

).5 

0.1017 

0.1021 

0.1021 

0.1021 

1.75 

0.1026 

0.1042 

0.1031 

0.1032 

1.0 

0.1036 

0.1022 

0.1042 

0.1043 


Parameter Ka 


2 „ 2 

(Ti , <72 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0101 

0.0103 

0.0103 

0.0103 

0.75 

0.0103 

0.0104 

0.0104 

0.0104 

1.0 

0.0103 

0.0102 

0.0105 

0.0106 


Parameter K3 


2^2 

,02 

0.5 

0.75 

1.0 

2.0 

0.5 

0.1019 

0.1023 

0.1022 

0.1023 

J.75 

0.1028 

0.1047 

0.1034 

0.1035 

1.0 

0.1039 

0.1025 

0.1046 

0.1047 


Parameter K4 


2 „ 2 
ai ,G2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0104 

0.0106 

0.0106 

0.0106 

0.75 

0.0107 

0.0112 

0.0109 

0.0109 

1.0 

0.011 

0.0106 

0.0112 

0.0112 


Results for Sample Size: 19 


Parameter Ki 


2 2 

1 » <^2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.1011 

0.1018 

0.1018" 

0.1018 

0.75 

0.1023 

0.1021 

0.1026 

0.1026 

1.0 

0.1031 

0.1018 

0.1028 

0.1034 


Parameter K3 


r 2 „ 2 

^1 >02 

0.5 

0.75 

1.0 

2.0 

0.5 

0.1012 

0.1019 

0.1019 

0.1019 

0.75 

0.1026 

0.1022 

0.1028 

0.1028" 

1.0 

0.1035 

0.102 

0.1029 

0.1036 


Parameter K2 


^ 2 „ 2 

CTl , CT2 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0101 

0.0102 

0.0102 

0.0102 

0.75 

0.0103 

0.0102 

0.0103 

0.0103 

1.0 

0.0104 

0.0102 

0.0103 

0.0103 


Parameter K4 


2 «. 2 
01 , 02 

0.5 

0.75 

1.0 

2.0 

0.5 

0.0103 

0.0105 

0.0105 

0.0105 

0.75 

0.0106 

0.0106 

0.0108 

0.0108 

1.0 

0.0108 

0.0105 

0.0108 

0.011 


A(\ 
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Table 4.12 Coverage Probabilities of Asymptotic Confidence Intervals 


ults for Sample Size: 10 


Parameter Ki 


^ 2 

0.5 

0.75 

1.0 

2.0 


919 

940 

959 

971 

5 

886 

910 

910 

960 

) 

879 

906 

901 

940 


Parameter K 3 


r«^ 

b 

0.5 

0.75 

1.0 

2.0 

5 

933 

935 

964 

972 

'5 

895 

922 

908 

960 

0 

878 

908 

900 

944 


Parameter K 2 


Cf2^ 

0.5 

0.75 

1.0 

2.0 

0.5 

925 

937 

948 

952 

0.75 

899 

919 

911 

948 

1.0 

898 

920 

914 

930 


Parameter K 4 


CTi ,a2 

0.5 

0.75 

1.0 

2.0 

0.5 

914 

940 

958 

966 

0.75 

894 

908 

923 

958 

1.0 

881 

904 

912 

938 


suits for Sample Size: 15 


Parameter Ki 


<I2^ 

0.5 

0.75 

1.0 

2.0 

5 

951 

962 

975 

995 

75 

925 

939 

947 

980 

.0 

917 

922 

948 

970 


Parameter K 2 


^ 2^2 
ai , cy2 

0.5 

0.75 

1.0 

2.0 

0.5 

951 

951 

948 

975 

0.75 

938 

935 

955 

955 

1.0 

934 

937 

95'0 

962 


Parameter K 3 


02 ^ 

0.5 

0.75 

1.0 

2.0 

J5 

948 

960 

972 

995 

75 

924 

938 

951 

979 

.0 

912 

921 

942 

964 


Parameter K 4 


^ 2 „ 2 

CTi , 02 

0.5 

0.75 

1.0 

2.0 

0.5 

953 

960 

973 

997 

0.75 

924 

943 

950 

984 

1.0 

917 

924 

948 

974 


^ults for Sample Size: 19 


Parameter Ki 


,CT2^ 

0.5 

0.75 

1.0 

2.0 

1.5 

951 

960 

964 

985 

.75 

929 

941 

955 

977 

1.0 

909 

927 

936 

975 


Parameter K 3 


• «• 2 
,cr2 

0.5 

0.75 

1.0 

2.0 

9.5 

950 

964 

965 

985 

1.75 

927 

942 

956 

980 

1.0 

907 

923 

940 

971 


Parameter K 2 


2 „ 2 
01 , 02 

0.5 

0.75 

1.0 

2.0 

0.5 

947 

952 

958 

952 

0.75 

938 

957 

953 

961 

1.0 

926 

938 

930 

955 


Parameter K 4 


^2^2 
cfi , cr2 

0.5 

0.75 

1.0 

2.0 

0.5 

951 

962 

970 

990 

0.75 

921 

945 

962 

977 

1.0 

914 

918 

935 

978 
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Table 4.13 Coverage Probabilities of Bootstrap Confidence Intervals 


Results for Sample Size: 10 


Parameter Ki 


CTi , CT 2 

0.5 

0.75 

1.0 

2.0 

0.5 

933 

926 

943 

934 

0.75 

939 

933 

932 

933 

1.0 

927 

943 

934 

934 


Parameter K3 


2 2 

ai , CT 2 

0.5 

0.75 

1.0 

2.0 

0.5 

941 

928 

943 

934 

0.75 

946 

942 

928 

939 

1.0 

925 

944 

934 

941 


Results for Sample Size: 15 


Parameter Ki 


^ 2 „ 2 
ai , a 2 

0.5 

0.75 

1.0 

2.0 

0.5 

937 

950 

936 

950 

0.75 

942 

940 

929 

943 

1.0 

950 

937 

949 

938 


Parameter K3 


2 „ 2 
ai ,02 

0.5 

0.75 

1.0 

2.0 

0.5 

938 

948 

937 

949 

0.75 

941 

940 

935 

941 

1.0 

948 

937 

949 

940 


Results for Sample Size: 19 


Parameter Ki 


^ 2^2 
CTl , CT 2 

0.5 

0.75 

1.0 

2.0 

0.5 

936 

926 

925 

936 

0.75 

924 

932 

920 

925 

1.0 

926 

923 

934 

935 


Parameter K2 


2 2 

Cl ,02 

0.5 

0.75 

1.0 

2.0 

0.5 

933 

925 

932 

942 

0.75 

931 

933 

928 

943 

1.0 

927 

942 

944 

924 


Parameter K4 


.-r 2 „ 2 
Cl , O2 

0.5 

0.75 

1.0 

2.0 

0.5 

908 

904 

919 

926 

0.75 

906 

910 

920 

919 

1.0 

906 

921 

927 

908 


Parameter K2 


^2^2 
Cl ,C2 

0.5 

0.75 

1.0 

2.0 

0.5 

931 

932 

940 

945 

0.75 

940 

929 

941 

935 

1.0 

939 

936 

951 

934 


Parameter K4 


... 2 _ 2 
Cl ,02 

0.5 

0.75 

1.0 

2.0 

0.5 

955 

964 

955 

959 

0.75 

956 

947 

951 

957 

1.0 

962 

954 

961 

956 


Parameter K2 


Cl ,C2 

0.5 

0.75 

1.0 

2.0 

0.5 

929 

914 

940 

909 

0.75 

924 

938 

927 

927 

1.0 

919 

939 

924 

928 


Parameter K3 


«■ 2 „ 2 

Cl , C 2 

0.5 

0.75 

1.0 

2.0 

0.5 

933 

925 

929 

941 

0.75 

922 

920 

920 

921 

1.0 

925 

929 

940 

933 


Parameter K4 


Cl ,C2 

0.5 

0.75 

1.0 

2.0 

0.5 

935 

932 

930 

941 

0.75 

917 

935 

932 

919 

1.0 

931 

932 

942 

933 
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Quantiles of Parameter K1 Sa 



Fig 4.3 Histogram & Q-Q Plot for Parameter Ki (n=10) 
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Quantiles of Parameter K2 Sample 


Histogram for Parameter K2 
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QQ Plot of K2 Sample Data versus Standard Normal 
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Quantiles of Parameter K2 Sample 


RESULTS & DISCUSSIONS 


Histogram for Parameter K2 



Standard Normal Quantiles 

Fig 4.8 Histogram & Q-Q Plot for Parameter K 2 (n=15) 
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4.2.4 Model Selection 

To see whether the model, from which we generated the data turns out to be 
true, the task of model discrimination has been performed. The cross validation 
method is used for model discrimination. The estimates of parameters obtained by 
both the methods of parameter estimation have been used for model selection. In the 
cross validation method the model, which has the smallest prediction error is 
considered as the best one. The value of the cross validation index was calculated 
using the Equation 3.66 & 3.67. The results are shown in the following Table 4.5. 
Table 4.14 Cross Validation Results 


mam 

CV(1) 

CV(2) 

CV(3) 

CV(4) 

0.5, 0.5 

21.5970 

25.9055 

3554.9 

21.6166 

0.5, 0.75 

22.3985 

29.3677 

1414.2 

22.4418 

0.5, 1.0 

31.5822 

33.6906 

1463.3 

32.9154 

0.5, 2.0 

42.3259 

45.0432 

1501.6 

43.9629 


From the above results it is explicit that for different values of variances, the 
cross validation index for model 1 is the least among all the rival models. This shows 
that model 1 is the best among all the rival models. It is tme since we generated the 
data from model 1 with different variances. The studies made by Buzzi et. al.[10] also 
says that model 1 is the best one with the help of posterior probabilities obtained by 
sequential design of experiments. 

4.3 Application to the real example 

Carr [11] discussed the catalytic isomerization of n-pentane and gave the 
kinetics for the reaction. Using the data given, die parameters were estimated by 
taking it as a single response system. 

4.3.1 Single response system 

The reaction scheme, which is explained by Carr [11] for the isomerization of 
n-pentane, involves two mechanisms. They are single site mechanism and dual site 
mechanism, which will be discriminated with the help of parameter estimation. The 
reaction is given as follows. 

(4.3) 
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For single site mechanism (Model 1) the rate equation will be represented by, 

tA(C,-C,/1.63) 

(l + t,C,+t,C,+t,C3)' 

For dual site mechanism (Model 2) the rate equation is given by, 

' KK{C2-cj\.6i) 


(4.4) 


(4.5) 


In the present work the theory of least square estimation is used for the estimation of 
the parameters involved in each of the response. The Objective functions are defined 
as the residual sum of squares. Considering the same kinetics as given in the literature 
the two objective functions are defined as 


(O.F.), = MinL 

.^3 


k,k,{c,-cj\.6i) 

^ (l + A^jC| + ^2^*2 + ^3^3 ) 



5 


(0-F\ = MinL 

*«.*1.*2.*3 


k,k,{c,-cJ\.e^) 


Y 

-Td ■ 


(4.6) 


(4.7) 


4.3.2 Parameter Estimation 

The best values of the parameters were calculated using a non-linear 
regression algorithm based upon the fimction iriiiiiniization technique of GNLM. The 
parameters obtained are compared with the values given by Jhonson et. al.[29] as 
follows . 


Table 4.15 Comparison of the Estimated Parameters 


Parameters 

Model 1 

Model 2 1 

Present 

Study 

Jhonson 
et. al.l291 

Present 

Study 

Jhonson 

et.al.r291 

Ao 

35.9071 

36.21 

133.0124 

132.97 

h 

0.0709 

0.9 

0.0022 

0.033 

k2 

0.0378 

0.48 

0.0013 

0.019 

ks 

0.1671 

2.13 

0.0049 

0.072 

SOS 

3.235 

4.31 

3.478 

4.88 


The variance of the above models is calculated by using. 


^^2 

CT 


SOS 

d.f. ■ 


(4.8) 
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QQ Plot of kO Sample Data versus Standard Normal 
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Quantiles of Parameter k3 Sample 
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Histogram for Parameter k3 



Standard Normal Quantiles 

Fig. 4.18 Histogram and Q-Q Plot for Parameter A 3 
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Standard Normal Quantiles 


Fig. 4.20 Histogram and Q-Q Plot for Parameter ki 
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-1 0 1 
Standard Normal Quantiles 


Fig. 4.22 Histogram and Q-Q Plot for Parameter *3 
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CHAPTER 5 

CONCLUSIONS AND RECOMMENDATIONS 

5.1 Conclusions 

In this dissertation some aspects of parameter estimation and model selection 
have been studied for single response and multi response systems, first with the 
simulated data and then with real life examples. The GNLM method is used for 
minimization of the residual sum of squares based on Asymptotic theory and 
Bootstrapping Methods. The program converges to minima that depend on the initial 
guesses, which should be taken properly in order to get the global minima. The case 
study of the single response system has been undertaken with a view to comparing the 
estimates of parameters and residual sum of squares with those reported in the 
literature. 

Based on this study the following conclusions can be drawn. 

1. In order to obtain the global minima, different initial guesses are used in the 
GNLM method and they are found to give same estimates. 

2. The efficiency of Least Square Estimates (LSE) is tested with the help of 
asymptotic theory and found that the assumptions of this theory are satisfied for 
large sample sizes (n > 15) . 

3. The average estimates and their mean squared errors give an insight into the 
robustness of the parameters. 

4. The asymptotic confidence intervals and coverage probabilities give an assurance 
to the 95% confidence interval of the parameters. 

5. More over bootstrapping confidence intervals and coverage probabilities works 
quite well even for lower sample sizes. 

6. From the normality plots of parameters obtained by simulation studies with 
different variances says the parameters are unbiased and normally distributed. 

7. Cross validation procedure can be used as a model selection technique in 
competing situations. 
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5.2 Recommendations 

In the minimization of the objective function, the initial guesses seem to be 
very important. Appropriate initial guesses lead to the identification of the global 
minima. The estimation of the parameters first with some simpler methods (for 
example by transferring the model to a linear one if possible) is recommended as 
these estimates of parameters can be used as the initial guesses for the methods in 
which the initial guesses are required. 

It is also recommended that one can use evolutionary algorithms like GA and 
SA, which does not require initial guesses but bounds of variables. The disadvantage 
of these methods is they are very time consuming. The combination of Asymptotic 
theory and Bootstrapping methods with Evolutionary algorithms can be tried to study 
the behavior of parameter estimation in nonlinear regression. 

Study of parametric and non-parametric bootstrapping methods on parameter 
estimation is recommended as the future work. 


70 



BIBLIOGRAPHY 


BIBLIOGRAPHY 


1. Bard Y. and Lapidus L., Kinetics analysis by digital parameter estimation, 
Catalysis Rev., 1968, 2, 67-73. 

2. Bard Y., Nonlinear Parameter Estimation, Academic Press, New York 1974. 

3. Bates D.M. and Watts D.G., Nonlinear Regression Analysis and Its Applications, 
Wiley, New York 1998. 

4. Bickel and Freedman, Asymptotic normality and the bootstrap in stratified 
sampling, Statist., 1984, 12,470-482. 

5. Boag I.F., Bacon D.W. and Downie J., Using a statistical Multi response method 
of experimental design in a reaction network study. Can. J. Chem. Eng., 1978, 56, 
389-395. 

6. Box G. E. P. and Draper N. R., The Bayesian estimation of common parameters 
from several responses, Biometrika, 1965, 52, 3, 355-365. 

7. Box G.E.P., Hxmter W.G. and Hunter J.S., Statisitics for Experimenters, John 
Wiley and Sons, 1978. 

8. Briemn, The little bootstrap and other methods for dimensionality selection in 
regression,/. Am. Statist. Assoc., 1992, 87,738-754. 

9. Buzzi F. G. and Forzatti P., A new sequential experimental design procedure for 
discriminating among rival models, Chem. Eng. Set, 1983, 38, 225-232. 

10. Buzzi F. G., Forzatti P., Emig G. and Hofinann H., Sequential experimental 
design for model discrimination in the case of multiple responses, Chem. Eng. 
Sci.,\9M, 39,81-85. 

11. Carr N.L., Kinetics of catalytic isomerization of n-Pentane, Ind. Eng. Chem., 
1960, 52, 391-399. 

12. Chemick M.R., Bootstrap Methods, A Practitioner’s guide, Willey series in Prob. 
&Stat.,l999. 

13. Cox D.R., Tests of separate families of hypotheses, Proc. Fourth Berkeley Sym., 
1961, 1, 105. 

14. Cox D.R., Further results on tests of separate families of hypotheses, J. Roy. Stat. 
Soc. , Series B, 1962, 24, 406. 


71 



BIBLIOGRAPHY 


15. Draper N. R. and Hunter G. W., Design of experiments for parameter estimation 
in multi response situations, Biometrika, 1961,53, 3, 525. 

16. a) Draper N. R. and Smith H., Applied Regression Analysis, John Wiley & Sons, 
Inc., Second Edition, 1981. 

b) Draper N. R. and Smith H., Applied Regression Analysis, John Wiley & Sons, 
Inc., Third Edition, 1998. 

17. Efron, The Jackknife, The Bootstrap and the other resampling plans, SIAM, 

1982a. 

18. Efron, Bootstrap confidence interval for a class of parametric problems, 
Biometrika, 1985, 72, 45-58. 

19. Efron, Better bootstrap confidence intervals, J. Am. Statist. Asoc., 1987, 82,171- 
200 . 

20. Efron B. and Tibishirani R., Bootstrap methods for standard errors: Confidence 
interval and measures of statistical accuracy. Statist. Sci., 1986, 1, 54-77. 

21. Efron B. and Tibishirani R., An Introduction to the Bootstrap, Chapman & Hall, 

1993. 

22. Freedman, Bootstrapping regression models, Statist., 1981, 9, 1218-1228. 

23. Freedman and Peters, Bootstrapping a regression equation: Some empirical 
results, J. Am. Statist. Assoc., 1984a, 79, 97-106. 

24. Froment G. F. and Bischoff K. B., Chemical Reactor Analysis and Design, John 
Wiley and Sons, Second edition, 1979. 

25. Hayes and Perl, Application of bootstrap statistical method to tau-decay-mode 
problem, Phys. Rev. D, 1989, 39, 274-279. 

26. Hill W. J., Statistical techniques for model building, PhD. Thesis, University of 
Wisconsin, 1966. 

27. Hocking R.R., The analysis and selection of variables in linear regression. 
Biometrics, 1976, 32, 1-49. 

28. Hunter W.G. and Reiner A.M., Designs for discriminating between two rival 
models. Technometrics, 1965, 7, 307. 

29. Jhonson R.A., Standal N.A. and Mezaki R., Re-Examination of Hougen-Watson 
rate model including a weighted least squares approach, I & EC Fundamentals, 
1968,7, 181-184. 


72 



BIBLIOGRAPHY 


30. Kanodia R., Parameter estimation and model discrimination in single response 
models, M Tech Thesis, I.I.T. Kanpur, 2001 . 

31. Kundu D. and Murali G., Model selection in linear regression, Comp. Statist. & 
Data Anal, 1995, 22, 461-469. 

32. Lehmann E.L. and Casella G., Theory of Point Estimation, Springer-Verlog, New 
York, Second Edition, 1998. 

33. Mezaki R. and Butt J. B., Estimation of rate constants from multi response kinetic 
data, Ind. Eng. Chem. Fund., 1968, 7 , 120. 

34. Moros R., Kalis H., Rex H.G. and Schaffarczyk St., A genetic algorithm for 
generating initial parameter estimations for kinetics of catalytic processes. Comp. 
& Chem. Engg, 1996, 20, 1257-1270. 

35. Negi A., Parameter estimation and model discrimination in single and multi 
response models, M. Tech. Thesis, LIT. Kanpur, 2003 . 

36. Prasad K.B.S. and Rao M.S., Use of expected likelihood in sequential model 
discrimination in multi response systems, Chem. Engg. Sc., 1976, 32, 141 1-1418. 

37. Rohatgi V.K., Statistical Inference, Wiley Series, New York, 1984. 

38. Seber G. A. F. and Wild, Non-Linear Regression Analysis, John Wiley and Sons, 
Inc., 1977. 

39. Seinfeld J.H., Nonlinear estimation theory, Ind. Engg. Chem., 1970, 62, 32-41. 

40. Sen and Srivastava M., Regression Analysis Theory, Methods and Applications, 
Springer-Verlag, New York, 1990. 

41. Shao J., On resampling methods for variance and bias estimation in linear models, 
Ann. Statist., 1988a, 16, 986-1008. 

42. Shao J., Bootstrap variance and bias estimation in linear models. Can. J. Statist., 
1988b, 16, 371-382. 

43. Shao J. and Tu. D., The Jackknife and Bootstrap, Springer-Verlag, New york, 
1995. 

44. Shimbakuro F.L, Lazer S., Dysen H.B. and Chemick M.R., A quasi-optical 
method for measuring the complex permittivity of materials, IEEE Trans., 
Microwave Theory Tech, 1984, 32, 661-665. 


73 


BIBLIOGRAPHY 


Ji5. Singh S., Design and analysis of experiments for model discrimination in 
uniresponse and multiresponse systems, Ph.D. Thesis, LIT. Kanpur, 1986. 

Stein K., Fitting kinetic equations to rate data assuming more than one set of 
active centers on a catalyst, I & EC Fundamentals, 1967, 6,2,169-174. 

,/n. Thompson M.L., Selection of variables in multiple regression, part I & II, Int. 
Statist. Rev., 1978, 1-19, 126-146. 

y4S. Vajda S. and Valko Y.P., A direct-indirect procedme for estimation of kinetic 
parameters. Comp. & Chem. Engg., 1986, 10, 49-58. 

. 49. Weber, On resampling techniques for regression models. Statist. Probab. Lett., 
1984, 2, 275-278. 

/50. Wu, Jackknife, Bootstrap and Other resampling plans in regression analysis, Ann. 
Statist., 1986, 14, 1261-1350. 



Response (Y2) Response (YQ 


APPENDIX I 


The simulated multi response data for different variances in the absence of co- 
variance is given in the following table. 



1 

12.5000 

11.9975 

11.7561 

12.1388 

12.8741 

13.2413 

13.6406 

12.1076 

2 

14.2857 

13.4588 

14.5418 

14.4518 

14.5286 

13.5488 

14.4832 

15.4977 

Dl 



18.1558 

17.7591 

18.7179 

18.8849 

16.9437 

17.7817 

4 

20.9302 

20.4491 

20.0622 

20.9449 

20.4410 

21.5367 

20.8006 

20.9586 

5 

16.6667 

15.4464 

16.4721 

16.9633 

17.1476 

17.3248 

16.2294 

17.3584 

6 

10.2740 

10.8490 

10.1605 

11.1162 

9.5156 

8.7905 

8.7765 

11.0506 

7 

22.7273 

23.7468 

21.9611 

23.2726 

23.3635 

24.0584 

22.3657 

22.8467 

8 

iH9ll 

EEmai 

12.2545 

11.7666 

12.1351 

12.2422 

13.1038 

12.2341 

9 

18.4211 

18.5191 

17.7779 

18.6229 

18.6224 

19.3819 

19.0453 

16.9951 

10 

26.4206 

25.8128 

26.4167 

27.0048 



26.7145 

26.7919 

11 

42.9078 

42.3759 

41.6891 

42.9021 

42.3609 

43.2416 

43.6724 

42.2693 

12 

22.3372 

23.2066 

23.2303 

22.9442 

22.4445 

23.1731 

22.5833 

22.8264 

mm 

27.9365 

28.7502 

27.5119 

28.4844 

27.6983 

28.6452 

27.8105 

26.7378 

14 

1.6129 

1.1830 

0.1535 

2.5363 

2.2993 

1.6450 

1.5353 

2.0115 

15 

13.4146 

13.9847 

13.4931 

14.2854 

13.3388 

12.6127 

13.6651 

13.9597 

16 

4.1985 

4.3520 

5.2504 

4.8763 

4.9151 

4.1015 

5.2072 

4.5274 

mm 

42.9078 

42.6437 

42.9453 

41.7379 

42.5717 

42.2124 

43.6056 

42.0302 

HB 

22.3372 

21.7488 

22.4517 

21.6368 

22.3859 

20.8952 

21.1631 

22.0364 

19 

27.9365 

28.1394 

27.9175 

28.4210 

28.2184 

27.1559 

28.0412 

28.3192 

1 

1.2500 

- 0.0362 

1.4003 

0.2751 

2.8287 

1.6660 

1.0576 

- 0.3402 

2 

1.4286 

0.3163 

2.1926 

0.8218 

2.3060 

1.8204 

1.9237 

0.7272 

3 

1.8182 

3.2435 

1.9761 

2.5050 

14113 

1.5246 

1.5583 

0.7406 

4 

2.0930 

2.0421 

2.7471 

2.1131 

0.1530 

2.1367 

3.0753 

3.0952 

5 

1.6667 

3.5256 

2.1066 

2.7305 

0.6967 

1.9901 

1.5113 

3.3961 

6 

1.0274 

0.8553 

1.1416 

- 0.3137 

1.4965 

1.1681 

- 0.2431 

1.7364 

7 

2.2727 

2.3953 

2.5153 



2.4548 

3.4811 

15248 

8 

1.3636 

2.0164 

0.5125 

- 0.2703 

1.7758 

2.8349 

1.7454 

1.5925 

9 

1.8421 

17158 

1.0245 

0.3994 

3.4078 

0.2319 

2.3317 

1.6186 

10 

2.6421 

2.2732 

2.6308 

2.9358 

2.9885 

2.8818 

2.0414 

1.7888 

11 

4.2908 

5.3033 

4.5977 

4.1504 

4.5241 

4.4958 

5.0129 

4.6364 

12 

2.2337 

1.6185 

1.4589 

1.1034 

2.8082 

2.7020 

0.2961 

2.3435 

13 

2.7937 

3.3647 

3.4970 

2.5011 

4.5133 

2.3829 

3.7442 

1.6606 

14 

0.1613 

- 0.1998 

0.2562 

- 0.4212 

2.2097 

0.7890 

0.1599 

- 0.5218 

15 

1.3415 

1.8672 

3.7071 

0.4451 

- 0.1083 

1.4630 


1.0636 

16 

0.4198 

1.0194 

0.7759 

0.6684 

0.7104 

1.0201 


1.0746 

17 

4.2908 

3.7040 

3.1590 

2.8011 

5.1236 

4.9723 

4.4704 

3.0424 

18 

2.2337 

2.6106 

2.5661 

2.5472 

1.8603 

3.1685 

2.4251 

1.6362 

19 

2.7937 

3.5240 

3.2262 

mm 

IIBl 

2.7482 

1.9224 

2.3118 


13.3516 

15.8118 

19.4180 

21.7198 

16.9497 

10.3343 

21.4284 

13.2742 

18.4028 

26.6184 

42.0347 

21.7616 

28.4199 

0.5836 

12.7430 

4.4332 

44.2371 

21.4259 

28.4783 

0.1220 

0.9852 

0.9665 

3.8736 

2.8807 

-19499 

1.7623 

2.1465 

- 0.3590 

2.3498 

3.6889 

2.9320 

1.5620 

0.2742 

0.6038 

-15797 

3.7473 

1.5861 

2.3814 
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Response (Y 2 ) Response (YQ 


S.No. 


iSSIPIl 

ESSS9BSQ 

1 

15.0807 

12.8273 

13.0551 13.1436 

2 

12.9795 

14.7614 

12.3281 15.3052 

3 

19.2053 

18.5806 

17.4214 19.1163 

4 

21.7080 

20.8575 

18.4863 22.1588 

5 

15.8328 

17.9815 

16.0073 16.4172 

6 

9.6873 

11.2523 

10.1592 9.5664 

7 

22.7929 

24.4494 

23.0274 22.1334 

8 

13.6240 

13.2241 

13.0524 13.3741 

9 

18.3441 

18.9862 

15.3473 19.6639 

10 

24.8620 

27.1605 

27.9717 24.8717 

11 

44.6104 

43.1279 

42.5004 42.5210 

12 

21.8682 

23.6500 

23.7653 22.6123 

13 

28.0311 

28.5657 

26.5833 28.7627 

14 

1.9000 

0.5049 

2.5169 0.6336 

15 

14.3340 

12.9676 

13.9563 13.3103 

16 

4.7086 

3.4725 

3.7335 4.3263 

17 

43.1532 

43.2618 

45.3382 42.9703 

18 

20.9366 

21.8304 

24.3576 22.7088 

19 

28.9061 

25.8328 

28.7338 27.8325 

1 

2.3769 

0.6744 

1.2810 0.2646 

2 

0.4118 

2.6844 

1.9693 0.8815 

3 

0.7333 

1.5326 

2.5021 1.8410 

4 

2.0402 

4.4329 

1.5030 4.0294 

5 

1.7243 

0.2508 

1.4056 2.2562 

6 

0.4311 

0.5629 

2.5445 1.1246 

7 

1.8735 

2.7466 

3.2800 2.6888 

8 

1.3437 

2.6563 

1.6671 2.0320 

9 

0.9627 

1.4479 

1.0250 4.3633 

10 

3.1604 

2.2122 

2.1509 3.0106 

11 

4.3329 

5.3605 

5.1582 6.4273 

12 

2.3392 

2.2692 

2.5945 2.6820 

13 

3.9221 

3.4418 

2.7133 3.9299 

14 

-0.3883 

1.2272 

0.9106 0.9828 

15 

2.4377 

1.6055 

-0.4505 3.8568 

16 

1.1658 

1.3532 

1.6331 -0.9068 

17 

4.1729 

4.9534 

ieb&eibb 

18 

2.4561 

1.0909 

ibseqiee^ 

19 

3.7974 

2.3527 

I 3.4032 I 4.0803 













































































































































































APPENDIX II 


The subroutine used for the minimization of objective function by GNLM 

method is given in this section. The code is written in MATLAB. 

function [beta,r,J] = gsnw (X, y, model, betaO) 
global s 

%function for Guass Newton-Levenberg Marquardt 
Method 

if (nargin<4) 

error ('GNLM requires four arguments.'); end 
if min{size(y)) ~= 1 

error ( ' Requires a vector second input 
argument . ' ) ; 
end 

y = y ( : ) ; 
a = size (y, 1) /2 ; 
si = sqrt (s (1) ) ; 
s2 = sqrt (s (2) ) ; 
bl = repmat (si, a, 1) ; 
b2 = repmat (s2, a, 1) ; 
b = [bl;b2]; 

if size(X,l) == 1 % turn a row vector into a column 
vector. 

X = X ( : ) ; 

end 

wasnan = (isnan(y) 1 any (isnan (X) , 2) ) ; 
if (any (wasnan) ) 
y (wasnan) = [ ] ; 

X (wasnan, : ) = [ ] ; 

end 

n = length (y); 
p = length (betaO) ; 
betaO = betaO ( : ) ; 

J = zeros (n,p) ; 

beta = betaO; 

betanew = beta +1; 

maxiter = 100; 

iter = 0; 

betatol = l.OE-4; 

rtol = l.OE-4; 

sse =1; 

sseold = sse; 

seps = sqrt(eps); 

zbeta = zeros (size (beta) ) ; 

slO = sqrt (10) ; 

eyep = eye(p); 

zerosp = zeros (p,l); 
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while (norm ( (betanew-beta) . / (beta+seps) ) > betatol i 
abs (sseold-sse) > rtol) & iter < maxiter 
if iter > 0, 

beta = be t anew; 

end 

iter = iter + 1; 

yfit = feval (model , beta, X) ; 

r = y - yfit; 

r = r . *b ; 

sseold = r ' *r; 

for k = l:p, 

delta = zbeta; 
if (beta(k) == 0) 

nb = sqrt (norm (beta) ) ; 
delta(k) = seps * (nb + (nb==0) ) ; 
else 

delta (k) = seps*beta ( k) ; 

end 

yplus = feval (model, beta+delta,X) ; 

J(:,k) = (yplus - yfit) /delta (k) ; 

end 

Jplus = [ J; (1. OE-2) *eyep] ; 
rplus = [r;zerosp]; 

% Levenberg-Marquardt type adjustment 
% Gauss-Newton step -> J\r 

% LM step -> inv ( J' *J+constant*eye (p) ) *J' *r 

step = Jplus\rplus; 

bet anew = beta + step; 

yfitnew = feval (model, betanew,X) ; 

rnew = y - yfitnew; 

rnew = rnew.*b; 

sse = rnew'*rnew; 

iterl = 0; 

while sse > sseold & iterl < 20 
step = step/slO; 
betanew = beta + step; 
yfitnew = feval (model, betanew, X) ; 
rnew = y - yfitnew; 
rnew = rnew.*b; 
sse = rnew' * rnew; 
iterl = iterl + 1; 

end 

end 

if iter == maxiter 

dispCGNLM did NOT converge. Returning results 
from last iteration.'); 
end 
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APPENDIX III 


The program used for the simulation of Single response system is given in this 
section. This program includes all the calculations of average estimates, mean 
squared errors, confidence intervals and coverage probabilities by asymptotic theory 
and bootstrapping methods. The code is written in MATLAB. 
clc; 

global k dt 
data; 

y = dt ( : , 3 ) ; 
c = dt ( : , 1 : 2 ) ; 
fp = fopen ( ' sim5 . txt ' , ' a ' ) ; 
fpl = f open ( ' dataS . txt ' , ' a ' ) ; 
sg = 1; 

%while sg<0.06 
xO = [400 5000] ; 
i = 0; 

pa(l:4) = 0; 
pb(l:4) = 0; 
while i <1000 
i = i+1 

er = randn(12,l); 
yt = y+er*sg; 

%Bootstrapping Perameters 
d = [yt c] ; 

bst = bootstrp(100,'m2fit',d); 
mbt(i,:) = mean (bst ); 

%Bootstrap Confidence Limits 

sbt = sort (bst); 

lb = (sbt (2, :)+sbt(3, :) )/2; 

ub = (sbt (97, : ) +sbt (98, : ) ) /2; 

cb = [ lb ' ub ' ] ; 

bt(i,:) = [cb(l,:) cb(2,:)]; 

k = 1; 

if bt (i,2*k-l)>400 | bt (i, 2*k) <400 
pb (k) = pb (k) +1; 

end 

k = k+1; 

if bt (i,2*k-l)>5000 | bt (i, 2*k) <5000 
pb(k) = pb(k)+l; 

end 

%Calculation of Perameters 
yt = yt ( : ) ; 

[x er J] = nlinf it (c, yt, @m2fun, xO) ; 
est (i,:)=x'; 

fprintf (fpl, ' %8 . 4f ', est(i,:)); 
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%Assymptotic Confidence Limits 

e = er'*er; 

a = diag ( inv ( J' * J) ) ; 

ns = size (yt, 1) ; 

b = 1 . 96 . *sqrt (a . *e/ (ns-3) ) ; 

cc = [x-b x+b] ; 

at(i,:) = [cc(l,:) cc(2,:)]; 

fprintf ( fpl , ' %8 . 4f at(i,:)); 

fprintf (fpl, ' \n ' ) ; 

k = 1; 

if at (i,2*k-l) >400 1 at (i, 2*k) <400 
pa (k) = pa (k) +1; 

end 

k = k+1; 

if at (i,2*k-l) >5000 | at (i, 2*k) <5000 
pa(k) = pa{k)+l; 

end 

end 

mt = mean (est ) ; 

fprintf (fp, ' \nMean Pearameters\n ' ) ; 
fprintf ( fp, ' %12 . 4f ',mt); 
mse = dia’g{cov(est)); 

fprintf (fp, ' \nMean Squared Errors of Perameters\n ' ) ; 
fprintf (fp, ' %12 . 4f ',mse); 
j = 0; 
while j < 2 

j = j+1; 

bi(:,j) = bt(;,2*j)-bt(:,2*j-l) ; 
ai(;,j) = at (: ,2*j)-at(: ,2*j-l) ; 
end 

mbi = mean (bi) ; 
mai = mean (ai) ; 

fprintf (fp, ' \nMean length of ACI\n'); 

fprintf (fp, ' %8 . 4f ',mai); 

fprintf ( fp, ' \nMean length of BtCI \n'); 

fprintf (fp, ' %8 . 4f ' ,mbi) ; 

pa = 1000-pa; 

fprintf (fp, ' \nCoverage Probabilities of ACI \n'); 
fprintf (fp, ' %4d ',pa); 
pb = 1000-pb; 

fprintf (fp, ' \nCoverage Probabilities of BtCI \n'); 
fprintf (fp, '%4d ',pb); 
k = 0; 
s = 100; 
qqpt (est (;,!)); 

str = [ ' qqp ' , int2str (s) , ' K' , int2str (k) ] ; 
saveas(gcf, str , 'fig') 
histf it (est ( ; , 1) , 30) 
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str = [ ' hst ' , int2str (s) , ' K' , int2str (k) ] ; 
saveas(gcf, str , 'fig') 
qqpt (est ( : , 2) ) ; 

str = [ ' qqp' , int2str (s) , ' K' , int2str (k) ] ; 
saveas(gcf, str , 'fig') 
histfit(est(:,2) ,30) 

str = [ ' hst ' , int2str (s) , ' K' , int2str (k) ] ; 
saveas(gcf, str , 'fig') 
fprintf (fpl, ' \n\n\n' ) ; 

%end 

fclose (fp) ; 

fclose (fpl) ; 

disp ( ' PROGRAMME END. ' ) ; 

function x=m2fit(dt) 
y = dt ( : , 1 ) ; 
c = [ dt ( : , 2 ) dt ( : , 3 ) ] ; 
xO = [400 5000] ; 

X = nlinfit (c, y, @m2fun, xO) ; 

X =x' ; 

function yhat=m2fun (s, c) 
n = ones (12,1) ; 

d = n+s (1) . *c( : , 1) . *exp(-s (2) . /c ( ; ,2) ) ; 
yhat = n./d; 
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APPENDIX IV 

The program used for the simulation of Multi response system is given in this 
section. This program includes all the calculations of average estimates, mean 
squared errors, confidence intervals and coverage probabilities by asymptotic theory 
and bootstrapping methods. The code is wntten in MATLAB. 

clc; 

global s k dt 
data; 

y = dt ( : , 3 : 4 ) ; 
c = dt ( : , 1 : 2 ) ; 
i i = 0.75; 

fp = fopen ( ' sim2 . txt ' , ' a' ) ; 
fpl = fopen { ' data3 . txt ' , ' a' ) ; 

jj = 1.0; 

%while jj < 2.1 
if jj == 1.25 
j j = 2.0; 

end 

s = [1/ii 1/j j ] ; 

mu = [0 0 ] ; 

cv = [ii 0; 0 j j ] ; 

fprintf (fp, ' \n\nFor variences: %4.2f %4 . 2f ' , ii, j j ) ; 
n = 10; 
while n<ll 
if n==20 
n=19; 

end 

fprintf (fp, ' \n\nResults for sample size: %d\n',n); 
cl = cn ( 1 : n, : ) ; 
y = yn (l:n, : ) ; 
c = repmat (cl, 2, 1) ; 
x 0 = [1 0.1 1 0 . 1 ] ; 
i = 0 ; 

pa (1:4) = 0; 
pb(l:4) = 0; 
while i <1000 
i = i+1; 

[n i] 

er = mvrnd (mu, cv, n) ; 
yt = y+er; 

%Bootstrapping Perameters 
d = [cl yt] ; 

bst = bootstrp(100,'crfit',d); 
mbt(i,:) = mean (bst); 

%Bootstrap Confidence Limits 
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sbt = sort(bst); 
lb = (sbt (2, : ) +sbt (3, : ) ) /2; 
ub = (sbt (97, ; ) +sbt (98, : ) ) /2; 
cb = [ lb ' ub ' ] /■ 

bt ( i , : ) = [ cb ( 1 , : ) cb ( 2 , : ) cb ( 3 , : ) cb ( 4 , : ) ] ; 
k = 1; 
while k<5 

if bt (i,2*k-l)>0.1 I bt (i, 2*k) <0. 1 
pb(k) = pb(k)+l; 

end 

k = k+1; 

if bt (i,2*k-l)>0.01 1 bt (i, 2*k) <0 . 01 
pb(k) = pb(k)+l; 

end 

k = k+1; 
end 

%Calculation of Perameters 
yt = yt ( : ) ; 

[x er J] = gsnw (c,yt, @plfun, xO) ; 
est (i, : ) =x'; 

fprintf (fpl, ' %8 . 4f ' , est(i,:)); 
f print f (fpl, ' \n' ) ; 

%Assymptotic Confidence Limits 
yl = plfun (x, c) ; 
el = yl-yt; 
e = el'*el; 
a = diag (inv ( J' * J) ) ; 
ns = size (yl, 1) ; 
b = 1. 96 . *sqrt (a. *e/ (ns-5) ) ; 
cc = [x-b x+b] ; 

at(i,:) = [cc(l,:) cc(2,:) cc(3,:) cc(4,:)]; 
k = 1; 
while k<5 

if at (i,2*k-l)>0.1 | at (i, 2*k) <0 . 1 
pa (k) = pa (k) +1; 

end 

k = k+1; 

if at (i,2*k-l)>0.01 | at (i, 2*k) <0 . 01 
pa(k) = pa(k)+l; 

end 

k = k+1; 
end 
end 

mt = mean (est); 

fprintf (fp, ' \nMean Pearameters\n ' ) ; 
fprintf (fp, ' %8 . 4f ' ,mt) ; 
mse = diag (cov (est) ) *le4; 

fprintf (fp, ' \nMean Squared Errors of PerametersXn ' ) ; 
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fprintf (fp, ' %12 . 4f ' ,inse) ; 
j = 0; 
while j < 4 

j = j+l; 

bi ( : , j ) = bt ( : , 2* j ) -bt ( : , 2* j-1) ; 
ai ( : , j ) = at'( : , 2* j ) -at ( : , 2* j-1) ; 
end 

mbi = mean(bi); 
mai = mean (ai) ; 

fprintf (fp, ' \nMean length of ACI\n'); 

fprintf (fp, ' %8 . 4f ' ,mai) ; 

fprintf (fp, ' \nMean length of BtCI \n'); 

fprintf (fp, ' %8 . 4f ' ,inbi) ; 

pa = 1000-pa; 

fprintf (fp, ' \nCoverage Probabilities of ACI \n'); 
fprintf (fp, ' %4d' , pa) ; 
pb = 1000-pb; 

fprintf (fp, ' \nCoverage Probabilities of BtCI \n'); 
fprintf (fp, ' %4d' ,pb) ; 
k = 0 ; 

qqpt (est ( : , 1) ) ; 

str = [ 'qqp' ,int2str(n) , 'K' ,int2str (k) ] ; 
saveas(gcf, str , 'fig') 
hist fit (est ( : , 1) , 30) 

str = [ ' hst ' , int2str (n) , ' K' , int2str (k) ] ; 
saveas(gcf, str , 'fig') 
qqpt (est ( : , 2) ) ; 

str = [ ' qqp' , int2str (n) , ' K' , int2str (k) ] ; 
saveas(gcf, str , 'fig') 
hist fit (est ( : , 2) , 30) 

str = [ 'hst' ,int2str(n) , 'K' ,int2str (k) ] ; 
saveas(gcf, str , 'fig') 
qqpt (est ( : , 3) ) ; 

str = [ ' qqp' , int2str (n) , ' K' , int2str (k) ] ; 
saveas(gcf, str , 'fig') 
hist fit (est ( : , 3) , 30) 

str = [ ' hst ' , int2str (n) , ' K' , int2str (k) ] ; 
saveas(gcf, str , 'fig') 
qqpt (est ( : , 4) ) ; 

str = [ ' qqp ' , int2str (n) , ' K' , int2str (k) ] ; 
saveas(gcf, str , 'fig') 
histfit (est ( : , 4) , 30) 

str = [ 'hst ' , int2str (n) , 'K' , int2str (k) ] ; 
saveas(gcf, str , 'fig') 
n = n+5; 

fprintf (fpl, ' \n\n\n' ) ; 
end 

jj = jj+0.25; 
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%end 

fclose(fp) ; 

fcloselfpl) ; 

disp ( ' PROGRAMME END. ' ) ; 


function c r f i t ( d ) 

global I.f:' 
dt = d; 

y = [dt ( : , 3) dt(:,4)]; 

y = y { : ) ; 

c = [dt(:,l) dt(:,2); dt(:,l) dt(:,2)]; 

xO = [0.1 0.01 0.1 0.01]; 

X = gsnw (c, y, §plfun, xO) ; 

X = X ' ; 

function [yhat,d]=plfun(s,cl) 

a = size (cl , 1 ) /2; 

c = cl (1 :a, : } ; 

xl = c ( : , 1) ; 

x2 = c ( : , 2 ) ; 

o = ones (size (xl) ) ; 

d = o+s (3) . *xl+s (4) . *x2; 

yl = (s(l) .*xl.*x2) ./d; 

y2 * (s (2) . *xl . *x2) . /d; 

yhat = {yl;y2]; 
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APPENDIX V 


The subroutine usc(i tbr the generation of multivariate random number is 
given here. The code i.s written in MATLAB. 

functi on X mvr nd Itnu, covm, n) 

I r.SVHK!) Gonorate multivariate normal random 

variab]<‘.'; . 

I 

% R - MVWJD{MU,COVM,N) Generates a sample of size 
N 

% random vntriablGs from the multivariate normal 
% distribution. MU is the d-dimensional mean as a 
% column vector. COVM is the d x d covariance 
matrix . 

if det (covm) <=^0 

% Then it is not a valid covariance matrix, 
error ('The covariance matrix must be positive 
def init.o’ ) 
end 

mu “ mu ( : ) ; % Just in case it is not a column 

vector. 

d “ length (mu); 

% get cholesky factorization of covariance 
R = chol(covm); 

% generate the standard normal random variables 
Z = randn (n, d) ; 

X * 2*R + ones (n, 1) *mu' ; 
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APPENDIX VI 


The program used for the cross validation procedure for single response 

system is given here. The code is written in MATLAB. 

clea r a 1 .1 

clc 

data 

dtl - lyl cl]; 
cv = 0 ; 

n=length (yl ) ; 
for i =- 1 : n 

ii = t (1: (i-1) ) C (i + 1) :n) ]; 
dt = dtl (ii, : ) ; 
y — dt ( • , 1 ) ; 
c ~ [dt ( : , 2 ) dt ( : , 3) ] ; 
xO - [400 5000] ; 

X - nlinf i t (c, y, @m3fun, xO) ; 

yhat “■ exp (“X ( 1) . ^cl (i, 1) . *exp(-x (2) . /cl (i,2) ) ) ; 
%yhat = 1/ (Hx{l 5 .*c:l (i,l) .*exp(-x(2) ,/cl(i,2) ) ) ; 
%yhat = (l+2*x(]).*cl{i,l).*exp(-x(2)./cl(i,2)))^ 

(-0.5) ; 

%yhat ( 1 « 3*x (1 ) . *cl (i, 1) . *exp (-X (2) . /cl (i, 2) ) ) '^ 
(-1/3) ; 

cv *= cv + (yhat-yl (i) ) ''2; 
i 

end 

cv 

function yhats=m3fun (s, c) 

Imodell 

yhat »= exp {-s { 1 ) . *c ( : , 1) • *exp (-S (2) . /c ( : , 2) ) ) ; 

%model3 & 4 

%n = ones (size (c, 1 ), 1) ; 

%d = nf 2 .*s(l) .*c(:,I) .*exp(-s(2) ./c(:,2)); 

lyhat « d. '' (-1/3) ; 
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APPENDIX VII 

The program used for the cross validation procedure for multi response 

system is given here. The code is written in MATLAB. 

global s 
clc 
data; 
ii = 0.5; 

jj = 2; 

n = 19; 

s = [1/ii 1/j j ] ; 

xO = [1 0.1 1 0.1] ; 

yt = yn; 

dtl = [yt cn] ; 

cv = 0; 

for i = l:n 

ii = [ (1: (i-1) ) ( (i+1) :n) ] ; 

dt = dtl (ii, : ) ; 
y = [dt(:,l) dt (:,2)] ; 
c = [dt ( : , 3) dt ( : , 4) ] ; 
c = r epmat ( c , 2 , 1 ) ; 
y = y ( : ) ; 

X = gsnw (c, y, @mrf 1, xO) ; 
dl = (1 + x(3) *cn(i,l)+x(4) *cn(i,2) ) ; 
d2 = (1 + X (3) *cn (i, 1) +x (4) *cn (i,2) ) ; 
yh(l) = (x(l)*cn(i,l)*cn(i,2)) /dl; 
yh(2) = {x(2)*cn(i,l)*cn(i,2))/d2; 
cv = cv + sum( (yh-yn(i, : ) ) .'^2) 
i 

end 

cv 

function yhat=mrf 1 (s, cl) 
a = size (cl , 1) /2 ; 
c = cl (l:a, : ) ; 
xl = c ( : , 1 ) ; 
x2 = c ( : , 2 ) ; 
o = ones (size (xl) ) ; 
d = o+s (3) . *xl+s (4) . *x2; 
yl = (s (1) .*xl.*x2) ./d; 
y2 = (s (2) .*xl.*x2) . /d; 
yhat = [yl;y2]; 
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